Graph Neural Preconditioners for Iterative Solutions of Sparse Linear Systems
Authors: Jie Chen
Venue: ICLR 2025 Poster (2025)
Keywords: General-purpose preconditioner
linear systems
graph neural networks
Primary area: other topics in machine learning (i.e., none of the above)
TL;DR: A new class of preconditioners based on graph neural networks is developed for reducing the condition number and speeding up the iterative solution of large sparse linear systems.

Abstract: Preconditioning is at the heart of iterative solutions of large, sparse linear systems of equations in scientific disciplines. Several algebraic approaches, which access no information beyond the matrix itself, are widely studied and used, but ill-conditioned matrices remain very challenging. We take a machine learning approach and propose using graph neural networks as a general-purpose preconditioner. They show attractive performance for many problems and can be used when the mainstream preconditioners perform poorly. Empirical evaluation on over 800 matrices suggests that the construction time of these graph neural preconditioners (GNPs) is more predictable and can be much shorter than that of other widely used ones, such as ILU and AMG, while the execution time is faster than using a Krylov method as the preconditioner, such as in inner-outer GMRES. GNPs have a strong potential for solving large-scale, challenging algebraic problems arising from not only partial differential equations, but also economics, statistics, graph, and optimization, to name a few.

Published as a conference paper at ICLR 2025
GRAPH NEURAL PRECONDITIONERS FOR ITERATIVE SOLUTIONS OF SPARSE LINEAR SYSTEMS
Jie Chen
MIT-IBM Watson Al Lab, IBM Research chenj jienjees jeuus. bmm com
ABSTRACT
Preconditioning is at the heart of iterative solutions of large, sparse linear sys- tems of equations in scientific disciplines. Several algebraic approaches, which access no information beyond the matrix itself, are widely studied and used, but ill-conditioned matrices remain very challenging. We take a machine learning approach and propose using graph neural networks as general-purposs precon- ditioner. They show attractive performance for many problems and can used when the mainstream preconditioners perform poorly. Empirical evaluation on over 800 matrices suggests that the construction time of these graph neural pre- conditioners (GNPs) is more predictable and can be much shorter than that of other widely used ones, such as ILU and AMG, while the execution time is faster than using Krylov method as the preconditioner, such as in inner-outer GM- RES. GNPs have a strong potential for solving large-scale, challenging algebraic problems arising from not only partial differential equations, but also economies, statistics, graph, and optimization, to name a few.
INTRODUCTION
Iterative methods are commonly used to solve large, sparse linear systems of equations Ax b. These methods typically build a Krylov subspace, onto which the original system is projected, such that an approximate solution within the subspace is extracted. For example, GMRES (Saad) & Schultz 1986), one of the most popularly used methods in practice, builds an orthonormal basis [v1, V2,.... J of the m-dimensional Krylov subspace by using the Arnoldi process and defines the approximate solution Xm Xo span((v1, Ym, Vm)), such that the residual norm l Axm//p 2 is minimized.
The effectiveness of Krylov methods heavily depends the conditioning of the matrix A. Hence, designing good preconditioner is crucial in practice. In some cases (e.g., solving partial differential equations, PDEs), the problem structure provides additional information that aids the development of a preconditioner applicable to this problem or a class of similar problems. For example, multigrid preconditioners are particularly effective for Poisson-like problems with mesh. In other cases, little information is known beyond the matrix itself. An example is the SuiteSparse matrix collec- tion (Davis& & Hul 2011). which contains thousands of sparse matrices for benchmarking numerical linear algebra algorithms. In these cases, general-purpose (also called "algebraic") preconditioner desirable; examples include ILU (Saad) 1994), approximate inverse (Chow & Saad 1998), and algebraic multigrid (AMG) (Ruge & Stüben 1987). However, it remains very challenging to design one that performs well universally.
In this work, we propose to use neural networks general-purpose preconditioner. Specifically, we consider preconditioned Krylov methods for solving problems in the form AMu b, where M s A- is the preconditioner, u is the new unknown, and Mu is the recovered solution. We design a graph neural network (GNN) (Zhou et al. 2020 Wu et al. (022) to assume the role of M and develop training method to learn M from select (b, pairs. This proposal is inspired by the universal approximation property of neural networks (Hornik et al. 1989 and encouraged by their
widespread success in artificial intelligence (Bommasani al. (2221).
Because a neural network is, in nature, nonlinear, our preconditioner M is not a linear operator anymore. We can longer build Krylov subspace for AM when a Krylov subspace is defined for

---

Published as a conference paper at ICLR 2025
only linear operators. Hence, we focus on flexible variants of the preconditioned Krylov methods instead. Specifically, we use flexible GMRES (FGMRES) (Saad [993). which considers M to be different in every Arnlli step. In this framework, all what matters is the subspace from which the
approximate solution is extracted. A nonlinear operator M can also be used to build this subspace and hence the neural preconditioner is applicable to FGMRES.
We use graph neural network because of the natural connection between a sparse matrix and the graph adjacency matrix, similar how AMG interprets the coefficient matrix with a graph. Many GNN architectures access the graph structure through A-mulipplications and they become polyno- mials of A when the nonlinearity is omitted (Wu al., [019 Chen al. 2020). This observation bridges the connection between a GNN as an operator and the polynomial approximation of A the latter of which an essential ingredient in the convergence theory of Krylov methods. In- terestingly, side benefit resulting from the access pattern of A is that the GNN is applicable to not only sparse matrices, but also structured matrices that admit fast A A-multiplications (such as hierarchical matrices (Hackbusch, 1999. Hackbusch & Börm 2002 Chen & Stein 2023). Toeplitz matrices (Chan & Jin 2007][Chen et al. (014). unassembled finite-element matrices (Ciarlet 2002).
and Gauss-Newton Newton matrices (Liu et al. 2023)).
There are a few advantages of the proposed graph neural preconditioner (GNP). It performs compar- atively much better for ill clldcondit problems, by learning the matrix inverse from data, mitigating the restrictive modeling of the sparsity pattern (as in ILU and approximate inverse), the insufficient quality of polynomial approximation (as in using Krylov methods the preconditioner), and the challenge of smoothing over a non-geometric mesh (as in AMG). Additionally, empirical evaluation suggests that the construction time of GNP is more predictable, inherited from the predictability of neural network training, than that of ILU and AMG; and the execution time of GNP is much shorter than GMRES as the preconditioner, which may be bottlenecked by the orthogonalization in the Amolli process.
1.1 CONTRIBUTIONS
Our work has a few technical contributions. First, we offer a convergence analysis for FGMRES. Although this method is well known and used in practice, little theoretical work was done, in part because the tools for analyzing Krylov subspaces are not readily applicable (as the subspace is no longer Krylov and we have no isomorphism with the space of polynomials). Instead, our analysis is posteriori, based on the computed subspace.
Second, we propose an effective approach to training the neural preconditioner; particular, train- ing data generation. While it is straightforward to train the neural network in an online-learning manner thrrough preparing streaming, random input- pairs in the form of (b, x) leaves open the question regarding what randomness best suits the purpose. We consider the bottom eigen- subspace of A when defining the sampling distribution and show the effectiveness of this approach.
Third, we develop scale-equivariant GNN as the preconditioner. Because of the way that training data are generated, the GNN inputs b can have vast different scales; meanwhile, the ground truth operator A 11 is equivariant the scaling of the inputs. Hence, to facilitate training, we design the GNN by imposing an inductive bias that obeys vcale-equivariance
Fourth, we adopt a novel evaluation protocol for general-purpol preconditioners. The common practice evaluates a preconditioner on a limited number of problems; in many occasions, these problems come from a certain type of PDEs or an application, which constitutes only a portion of the use cases. To broadly evaluate a new kind of purppss preconditioner and understand its strengths and weaknesses, we propose to test on large number of matrices commonly used by the community. To this end, we use the SuiteSparse collection and perform evaluation on a substantial portion of it (all non-SPD, real matrices falling inside a size interval; over 800 from 50 application areas in this work). To streamline the evaluation. define evaluation criteria, limit the tuning of hyperparameters, and study statistics over the distribution of matrices.
1.2 RELATED WORK
is important to put the proposed GNP in context. The idea of using a neural network to con- struct preconditioner emerged recently. The relevant approaches learn the nonzero entries of the
2

---

Published as a conference paper at ICLR 2025
incomplete factors (Li et al. 2023 Häusner et al. (023) or their correction Trifonov et al. 2024 or of the approximate inverse (Bänkestad et al. 2024). In these approaches, the neural network does not directly approximate the mapping from b to x like ours does; rather, it is used to complete the nonzeros of the predefined sparsity structure. A drawback of these approaches is that the predefined sparsity imposes nonzero lower bound the approximation of A A- which (and whose factors) are generally not sparse. Moreover, although the neural network can be trained on a PDE problem (by varying the coefficients, grid sizes, or boundary conditions) and generalizes well the same problem, it is unlikely that a single network can work well for all problems and matrices in life.
Approaches wherein the neural network directly approximates the matrix inverse are more akin to inysice- neural networks (PINNs) (Raissi et al. (2019) and neural operators (NOs) (Rudikov et al.) 2024] Kovachki et al., 2022 Li et al.J 2020; 2021] Lu et al., 2021). How- ever, the learning of PINNs requires a PDE whose residual forms a partial training loss, whereas NOs consider moinite-dimensionn
function spaces, with PDEs being the application. In our case of ceneral-purposs preconditioner, there is not an underlying PDE; and not every matrix problem relates to function spaces with additional properties to exploit (e.g., spatial coordinates, smoothness, and decay of the Green's 's function). For example, one may be interested in finding the commute times between a node and all other nodes in graph; this problem can be solved by solving a linear system with respect the graph Laplacian matrix. The only information exploit is the matrix itself.
2
METHOD
Let A E Rax", From now on, M is an operator I - IR" and we write M(v) to mean applying the operator on R". This notation includes the special case M(v) = Mv when M is matrix.
2.1 FLEXIBLE GMRES
A standard preconditioned Krylov solver solves the linear system AMu = b by viewing AM as the new matrix and building a Krylov subspace, from which an approximate solution Um extracted and Xm is recovered from Mu,. It is important to note that a subspace, by definition, is linear. However, applying a nonlinear operator M on the linear vector space does not result in vector space spanned by the mapped basis. Hence, convergence theory of xm
MM,,.) is broken.
We resort to flexible variants of Krylov solvers (aad 1993 Notay 2000 Chen et al. 2016) and focus on flexible GMRES for simplicity. FGMRES was designed such that the matrix M can change in each Arnoldi step: we extend it to a fixed, but nonlinear, operator M.
Algorithm Tin in Section B summarizes the solver. The algorithm assumes a restart length m. Starting from an initial guess Xo with residual vector rp and 2-norm B, FGMRES runs m Arnoldi steps, resulting in the relation
AZm Vn Hm + Am I.m m 1 V m H m:
(1)
where V.. Mi,.... vm] E Ruxm is orthonormal, Zm = [Z]..... Zm] E R Raxn contains m
columns z, M(v,), H. [byg] E mm is upper Hessenberg, and Hn, = [bal] E TR (m)+))x extends Hn, m. with an additional row at the bottom. The approximate solution x,, is computed in the form X., E Xg + span( (oaan(z)...... zm)) = {xo + Zn.)). J. such that the residual norm |b - Axm l|2 = ||ße| 112 is minimized. If the residual norm is sufficiently small, the solution is claimed; otherwise, let xm be the initial guess of the next Arnoldi cycle and repeat.
We provide a convergence analysis below. All proofs are in Section C
Theorem 1. Assume that FGMRES is run without restart and without breakdown. We use the rilde notation to denote the counterpart quantities when FGMRES is run with a fixed preconditioner ma- trix M. For any M such that AM can be diagonalized, as in X AMMX = = dia. Xy))
the residual r.m b Axm. satisfies
|mmll 1 S K2 (X)e (m)) (A))|roll2 1/Q... QT am Ô.h 11 |2llrol|2
where K2 denotes the 2-norm condition number, Qm (resp. Qm) is the thin- Q factor of AZ (resp. A AZm. Pm denotes the space of degree-m polynomials, and
(m) (A)= =
PEP .m: p()== Ip(Aa))
min max
3

---

Published as a conference paper at ICLR 2025
Let us consider the two terms on the right-hand side of (2) The first term N2 (X)) (A) |ro||| is the standard convergence result of GMRES on the matrix AM. The factor c(m) (A) ) is approximately exponential in m., when the eigenvalues X; are located in ellipse which excludes the origin; see [Sad (2003, Corollary 6.33). Furthermore, when the eigenvalues are closer to each other than to the origin, the base of the exponential is close to zero, resulting in rapid convergence. Meanwhile, when A MM is close to normal, K2g (X) is close to 1.
What M affects is the second term QQm Q.m. Q 12/.la//rolll This term does not show conver- gence on the surface, but one can always find an M such that it vanishes, assuming no breakdown. When M is close A A! AZn a small perturbation of V by definition. Then, suggests that Hm close to the identity matrix for all m and so is Hn.
Corollary 2. Assume that FGMRES is without restart and without breakdown. On completion, let H. be diagonalizable, as in Y Y = E = diag(o)......00 Then, the residual norm satisfies
|mmll S Kg Y) (m) (E)l|ollla
for all m.
See Section D for an approximately exponential convergence of [mmll due to ((m) (E).
Note that while our results assume infinite precision, another analysis (Arioli & Duff 2009) was conducted based on finite arithmetics. Such analysis uses the machine precision to bound the resid- ual norm for a sufficiently large m, asserting backward stability of FGMRES; the other hand, can obtain the convergence rate under infinite precision.
2.2 TRAINING DATA GENERATION
We generate set of (b, x) pairs to train M b - x. Unlike neural operators where the data generation is costly (e.g., requiring solving PDEs), which limits the training set size, for linear systems, we can sample x from a certain distribution and obtain b Ax with negligible cost, which creates streaming training set of (b, x) pairs without size limit
There is no free lunch. One may sample x ~ I,),,,) which leads to b AA J. This distribution is skewed toward the dominant eigen-subspace of AA Hence, using samples from it for training may result in a poor performance of the preconditioner, when applied to inputs lying close the bottom eigen-subspace. One may alternatively want the training data b ~ (O, I,), which covers uniformly all possibilities of the preconditioner input on a sphere. However, in this case, (O A -1A-T] is also skewed in the spectrum, causing difficulty in network training.
We resort to the Arnoldi process to obtain an approximation of A A which strikes a balance. Specifically,
the Amoldi process in m steps without a preconditioner, which results in simplification of []]
AV m Hm + hm Im V Vm- em
V V H
(4)
Let the singular value decomposition of Hm be W S,, ZT where W.. W. [(m+1)xm and Zm Rmxm are orthonormal and Sm E Rmxm is diagonal. We define
E
x V m Zm
E e~ I(,,
(5)
Then, with simple algebra we obtain that
x E range( mn) b E range( ( Wn+1),
x NO, )xx) b ~~~ N Eb)
X HV
Σ = (V (Vm))mm)) W m))
The covariance of b. Σb is projector, because V M and Wm m are orthonormal. As grows, this covariance tends to the identity I,, By the theory of the Arnoldi process (Morgan 2002). V Vm+1 contains better and better approximations of the extreme cigenvectors (0, Σb) has a high density close to the bottom eigen-subspace A. Meanwhile, because b is
increases. Hence,
Neural operators learn the mapping from the boundary condition to the solution. Strietly speaking, while is the solution, b is not the boundary condition.
4

---

Published as a conference paper at ICLR 2025
XL
MLP
GCONV
MLP
7X
xxd
nxd
nx1
Figure 1: Our GNN architecture for the preconditioner M. The nonlinearity is ReLU. The component denotes the scale-equivariant operator
low-rank, the distribution has zero density along many directions. Therefore, we sample x from both (O, ) and (0, I,) to form each training batch, to cover more important subspace directions.
note that the use of Gaussian distributions brings the convenience to generally analyze input/output distributions. It is possible that other distributions more effective if knowledge of exists and can be exploited.
2.3 SEALL-EQUIVARIAE GRAPH NEURAL PRECONDITIONER
We now define the neural network M that approximates admits graph interpretation, similar to how AMG interprets the coefficient matrix with graph. Treating as the graph adjacency matrix allows us to use GNNs (Kipf & Welling 2017] Hammlton al. 2017; Veličkovié et al.) 2018 Xu et al., 2019) parameterize the preconditioner. These neural networks perform graph convolutions in a neighborhood of each node, reducing the quadratic pairwise cost by ignoring nodes faraway from the neighborhood.
Naturally, a (sparse) matrix A [atj]
Our neural network based on the graph convolutional network (GCN) (Kipf & Welling, 2017). A
GCN layer is defined as GCONV(X) = ReLUU(AX), where Å E R is some normalization of A, X E Raxda is the input data, and W E R. Xdua is learnable parameter.
Our neural network, as illustrated in Figure I enhances from GCN in few manners:
1. We normalize differently. Specifically, we define
A = A/7 where min max, { E, |alss }. max,
Ex [als)
(6)
The factor y is an upper bound of the spectral radius of A based on the Gershgorin circle the- orem (Gerschgorin 1931). This avoids division-by-zero that may occur in the standard GCN normalization, because the nonzeros of A can be both positive and negative.
2. We add residual connection to GCONV, to allow stacking deeper GNN without suffering the smoothing problem (Chen al. 2020j. Specifically, we define Res- Res-GCONV layer to be
Res- Res-GCONV(X) ReLU(XU + AxW)
Moreover, we let the parameters U and W be square matrices of dimension d x d, such that the input/output dimensions across layers do not change. We stack L such layers.
We use an MLP encoder (resp. MLP decoder) before (resp. after) the L Res-GCONV layers to lift (resp. project) the dimensions.
Additionally, the main novelty of the architecture scale-equivarianc. Note that the operator to approximate, A-1, is scale-equivariant, but general neural network is not guaranteed so. Because the input may have very different scales (due to the sampling approach considered in the preceding subsection), to facilitate training, design a parameter-free scaling s and back-scaling
/ and --)
where |||l2
Vn
We apply at the beginning and s- at the end of the neural network. Effectively, restrict the input space of the neural network from the full R² to a more controllable space- the sphere - Clearly, guarantees scale-equivariance; i.e., M(ab) oM(b) for any scalar + O.

---

Published as a conference paper at ICLR 2025
3 EVALUATION METHODOLOGY AND EXPERIMENT SETTING
Problems. Our evaluation strives to cover large number of problems, which come from di- verse application areas as possible. To this end, we turn to the SuiteSparse matrix collection https://spares tamu edul, which is widely used benchmark in numerical linear algebra. We select all square, real-valued, and non-SPD matrices whose number of rows falls between 1K and 100K and whose number of nonzeros is fewer than 2M. This selection results in 867 matrices from 50 application areas. Some applications are involved with PDEs (such as computational fluid dynamics), while others come from graph, optimization, economics, and statisties problems.
To facilitate solving the linear systems, we prescale each A by y defined in (6). This scaling is a form of normalization, which avoids the overall entries of the matrix being exceedingly small or large. All experiments assume the ground truth solution x = 1 FGMRES starts with Xo O.
Compared methods. We compare GNP with three widely used general-purpose preconditioners ILU, AMG, and GMRES (using GMRES to precondition GMRES is also called inner-outer GMRES) and one simple-to-mmplemen but weak preconditioner- conditioning over 800 matrices from different domains limits the choices, but these preconditioners are handy with Python software support.
The endeavor of pre-
ILU comes from scipy. sparse. linalg.spilu; it in turn comes from SuperLU Li 2005
Li & Shao 2010), which implements the thresholded version of ILU (Saad, 1994). We use the default drop tolerance and fill factor without tuning. AMG comes from PyAMG (Bell jet al.] (023). Specifically, we use pyamg waamn.blackeen www.bbakoo............o.o....... sol solvere ver () .aspreoondit oner as the preconditioner, with the configuration of the blackbox solver computed from traag.blackboor
solver........ wwwver.configuree.........oo...e www.oooggggiiioo.ommom GMRES, we use 10 iterations and stop when it reaches a relative residual norm tolerance of 1e-6. For (the outer) FGMRES, the restart cycle is 10.
GMRES self-implemented. For (the inner)
Note that sophisticated preconditioners require tremendous efforts implement robustly. Hence, we prioritize ones that are more robust and that come with continual supports. For example, we also experiment with AmgX (Naumov et al. 2015) but find that while heing faster, it throws more erros than does PyAMG. See Section F.4 for details.
Stopping criteria. The stopping criteria of all linear system solutions rtol 1e-8 and maxiters = 100. For time comparisons, we additionally run experiments by enabling timeout and disabling maxiters. We the timeout to be the maximum solution time among all preconditioners when using maxiters = 100 as the stopping criterion.
Evaluation metrics. Evaluating the solution qualities and comparing different methods for large number of problems require programmable metrics, but defining these metrics is harder than it ap- pears to be. Comparing the number of iterations alone is insufficient, because the per-iteration time of each method is different. However, to compare time, one waits until the solution reaches rtol, which impossible to set if one wants all solutions to reach this tolerance. Thus, we can use maxiters and t timeout to terminate solution. However, one solution may reach maxi wiiteers early (in time) and attain poor accuracy, while another solution may reach timeout very late but attain good accuracy. Which one is better is debatable.
Hence, we propose two novel metrics, depending on the use of the stopping criteria. The first one, area under the relative residual norm curve with respect to iterations, is defined
iters
Iter-AUC-
log10 Ti
log rtol, T4 = lb Ax, |
1=0
where i iters is the actual number of iterations when FGMRES stops and r; is the relative residual norm iteration We take logarithm because residual plots are typically done in the log scale. This metric compares methods based on the iteration count, taking into account the history (e.g.. convergence speed at different stages) rather than the final error alone. This metric is used when the stopping criteria are rtol and maxiters.
TTii is common practice. Avoid setting to he Gaussian random, because this (partially) matches the training data generation, resulting the inability to test out-of-distribution cases.

---

Published as a conference paper at ICLR 2025
by #iterations
by solution time
by #iter
by time
25%
1.91c+00
1.18c+00
50%
6.74c+00 2. 33c+000
34.9%
31.5%
Best preconditioner
75%
6.78e+03 8. 36101
13.7%
19.4%
LU
1008
.89++10 1.91ee10
0.2%
3.7%
AMG
Table 1: Distribution of the
Jacobi
27.3%
2.1%
21.7%
27.6%
6.7%
11.2%
No precond
GMRES precond
residual-norm ratio between the second best preconditioner and GNP, when GNP performs the best. Distribution de- scribed by using percentiles.
Figure 2: Percentage of problems on which each preconditioner per- forms the best.
GNP
ILU
AMG
Jacobi
GMRES precand
Time (seconds)
10
10
dank
10-1
10-1
10
0
10 0
10
0
10
10
matrix 1ee+04
Construction
Solution
Figure 3: Preconditioner construction time and solution time (using maxiters to stop). The con- struction time of Jacobi is negligible and not shown. GMRES does require construction.
The second metric, area under the relative residual norm curve with respect time, is defined as
Time-AUC = [f [10810 r(1) - log10 rtol] dl Σ
iters
1ogi0 r; - log10 rtol][4; - 4411,;
i=1
where 7' is the elapsed time when FGMRES stops, r(t) is the relative residual norm at time t, and t; the elapsed time iteration This metric compares methods based on the solution time, taking into account the history of the errors. It is used when the stopping criteria are rtol and timeout.
Hyperparameters. To feasibly train over 800 GNNs, we use the same set of hyperparameters for each of them without tuning. There is strong potential that the GNN performance can be greatly improved with careful tuning. Our purpose is to understand the general behavior of GNP, partieu- larly its strengths and weaknesses. We use L 8 Ree--CONV layers, set the layer input/output dimension to 16, and use 2-layer MLPs with hidden dimension 32 for lifting/projection. We use Adam (Kingma & Ba 2015) as the optimizer, set the learning rate to 1 1e-3, and train for 2000 steps with a batch size of 16. We apply neither dropouts nor weight decays. Because data are sampled from the same distribution, resorting to a validation set.
use the model at the best training epoch as the preconditioner, without
We use the & residual norm AAM(b) - Ax|li the training loss, common practice of robust regression. We use m = 40 Arnoldi steps when sampling the (x, b) pairs according to (5) Among the 16 pairs batch, 8 pairs follow (5) and 8 pairs follow I.)
Compute environment. Our experiments are conducted on a machine with one Tesla V100(16GB) GPU, 96 Intel Xeon 2.40GHz cores, and 386GB main memory. All code is implemented in Python with Pytorch. All computations (training and inference) use the GPU whenever supported.
4 RESULTS
We organize the experiment results by key questions of interest.
How does GNP perform compared with traditional preconditioners? We consider three factors: convergence speed, runtime speed, and preconditioner construction cost.
In Figure 2 we show the percentage of problems on which each preconditioner performs the best, with respect to iteration counts and solution time (see the Iter-AUC and Time-AUC metrics defined
7

---

Published as a conference paper at ICLR 2025
120
precondiooncr
100
400
100
GNP
count
300
50
200
40
Jacobl
GMRES as pracond
20
No precond
10
Chem
matrix size/
logl0icondition numberr
Figure 4: Breakdown of best preconditioners with respect to matrix sizes, condition numbers, and application areas. Only the application areas with the top number of problems are shown. The last bar in the middle plot is for condition number 2 10
in the preceding section). GNP performs the best for a substantial portion of the problems in both metrics. In comparison, ILU and AMG perform the best for more problems, while Jacobi, GMRES (as preconditioner), and no preconditioner performs the best for fewer problems in the time metric.
That ILU and AMG are the best for more problems does not compromise the competitiveness of GNP. For one reason, they are less robust and sometimes bare a significantly higher cost construc- tion, will be revealed later. Another reason is that GNP can be used when other preconditioners are less satisfactory. In Table[] we summarize the distribution of the residual-neg norm ratio between the second best preconditioner and GNP, when GNP is the best. The ratio means factor of decrease in the residual norm if GNP is used in place of other preconditioners. Under the Iter-AUC metric, factor of 6780x decrease can be seen the 75th percentile, and in the best occasion, more than ten orders of magnitude decrease is seen. Similarly under the Time-AUC metric. These significant reductions manifest the usefulness of GNP.
In Figure 3 and Figure [ of Section F22 we plot the preconditioner construction time and solution time for each matrix and each preconditioner. Here, the solution is terminated by maxiters rather than t iimeuut (such that the times are different across matrices). Two observations can be made. First, the training time of GNP is nearly proportional the matrix size, while the construction time of ILU and AMG is hardly predictable. Even though for many of the problems, the time is only fraction of that of GNP, there exist quite a few cases where the time significantly longer. In the worst case, the construction of the AMG preconditioner is more than an order of magnitude more costly than that of GNP.
Second, there is clear gap between the construction and solution time for GNP, but such gap exists for ILU and AMG. One is tempted to compare the overall time, but one preconditioner can be used for any number of right-hand sides, and hence different conclusions may be reached regarding time depending on this number. When it is large, the solution time dominates, in which case one sees that most of the cases GNP is faster ILU, AMG, and GMRES.
On what problems does GNP perform the best? We expand the granularity of Figure into Figure overlaid with the distributions of the matrices regarding the size, the condition number, and the application area, respectively. One finding is that GNP is particularly useful for ill-conditioned matrices (i.e., those with condition number > 1016). ILU less competitive in these problems, possibly because it focuses on the sparsity structure, whose connection with the spectrum is less clear. Another finding that GNP is particularly useful three application areas: chemical process simulation problems, economic problems, and eigenvalue/moodl reduction problems. These areas constitute many problems the dataset; both the number and the proportion of them in which GNP performs the best are substantially high.
On the other hand, the matrix size does not affect much the GNP performance. One sees this more clearly when consulting Figure 9 Section F3 which is "proportion" version of Figure the proportion problems on which GNP performs the best is relatively flat across matrix sizes. Similarly. symmetry does not appear to play a distinguishing role either: GNP performs the best 7.6% of the symmetric matrices and 17.5% of the nonsymmetric matrices.
How robust is GNP? We summarize in Table the number of failures for each preconditioner. GNP and GMRES are highly robust. We see that neural network training does not cause any troubles, which is practical advantage. In contrast, II.U fails for nearly half of the problems, AMG fails for
8

---

Published as a conference paper at ICLR 2025
Table 2: Failures of preconditioners (count and proportion).
GNP
II.UU
AMG
Jacohi
GMRES as precond
Construction failure Solution failure
(0.00%)
348 (40. 14%)
62 (7.15%)
N/A
N/A 2 (0.23%)
(0.)
61 ( 7.04%)
5 (0.58%)
(6.11%)
G-15, indeffbratudd 2D/3D problemn cond 5.92e+02
VenVellee/std1 Jac3 chernical process simulation prablem cond 1.02e+22
Simmn/xenkat75 subsequent cfd problem cond NJA
Hollingerterrk3jacaa economic problem
cond 6.79e+
logibirelative residucl)
-3-
-3
-8
0
100
0
100
0
100
0
50
100
Iterations
GNP
AMG
jacobl
GMRES as precond
No precond
Figure 5: Example convergence histories.
nearly 8%, and Jacobi fails for over 6%. According to the error log, the common failures of ILU are that "(f)actor is exactly singular" and that "matrix is singular in file ilu dluddivvo...." while the common failures of AMG are "array ... contain(s) infs or NaNs". Meanwhile, solution failures occur when the residual norm tracked by the QR factorization of the upper Hessenberg H., fails match the actual residual norm. While ILU and AMG perform the best for more problems than does GNP, safeguarding their robustness is challenging, rendering GNP more attractive.
What does the convergence history look like? In Figure 5 we show a few examples when GNP performs the best. Note that in this case, ILU fails for most of the problems, in either the construction or the solution phase. The examples indicate that the convergence of GNP either tracks that of other preconditioners with a similar rate (S mimmoo//e venkat or becomes significantly faster. Notably, on Vanvel zen/ //sst_d_sas _Jac3, GNP uses only four iterations to reach 1,-8, while other preconditioners take 100 iterations but still cannot decrease the relative residual to below 1e-2.
We additionally plot the convergence curves with respect to time in Figure 10 in Section F.6 An important finding to note is that GMRES as a preconditioner generally takes much longer time to run than GNP. This is because GMRES is bottlenecked by the orthogonalization process, shadowing the trickier cost comparison between more matrix-vector multiplications in GMRES and fewer matrix- matrix multiplications in GNP.
Are the proposed training-data generation and the scale-equivariance design necessary? In Figure [6] we compare the proposed designs versus alternatives, for both the training loss and the relative residual norm. In the comparison of training data generation, using x ~ (o, I.) leads to lower losses, while using b (0, I,) leads to higher losses. This is expected, because the former expects similar outputs from the neural network, which is easier to train, while the latter expects drastically different outputs, which make the network difficult to train. Using x from the proposed mixture leads to losses lying in between. This case leads to the best preconditioning performance more problems having lower residual norms (especially those near rtol). other words, the proposed training data generation strikes the best balance between easiness of neural network training and goodness of the resulting preconditioner.
In the comparison of using the scaling operator versus not, we see that the training behaviors barely differ between the two choices, but using scaling admits an advantage that more problems result in lower residual norms. These observations corroborate the choices make in the neural network design and training data generation.
DISCUSSIONS AND CONCLUSIONS
We have presented a GNN approach to preconditioning Krylov solvers. This is the first work in our knowledge to neural network as an approximation the matrix inverse, without exploiting
9

---

Published as a conference paper at ICLR 2025
0.24
1.50
0.21
140
0.20
Vatient scaln-egutvarte
Density
1.20
sormal
0.10
0.15
normal
100
1.20
0.16
0.90
0.12
0.80
0.12
0.60
0.09
0.60
0.08
-20
-10
-2.5
0.0
2.5
-20
-10
log10ltraining less)
loglo(relative residual)
loq10ltraining loss)
loglO(relative residual)
Figure 6: Left: comparison of training data generation; right: comparison of scale-equivariance.
information beyond the matrix itself (i.e., an algebraic preconditioenr). Compared with traditional algebraic preconditioners, this approach is robust with predictable construction costs and is the most effective many problems. We evaluate the approach on more than 800 matrices from 50 applica- tion areas, the widest coverage in the literature.
One common skepticism on neural network approaches for preconditioning is that in general, net- work training is costly and in our case, the network can be used for only one A. We have shown to the contrary, our training cost sometimes can be much lower than the construction cost widely used preconditioners such as ILU and AMG, because network training is more predictable while the other preconditioners suffer the manipulation of the irregular nonzero patterns. Moreover, a benefit of neural networks is that compute infrastructures and library supports are in place, which facilitate implementation and offer robustness, as opposed to other algebraic preconditioners whose imple- mentations challenging and robustness is extremely difficult to achieve. Additionally, traditional algebraic preconditioners are constructed for only one A as well, not more flexible than ours.
While one may hope that a single neural network can solve many, if not all, linear systems (different A's), as recent approaches (PINN, NO, or learning the incomplete factors) appear to suggest, it is unrealistic expect that a single network has the capacity to learn the matrix inverse for all matrices. Existing approaches work on a distribution of linear systems for the same problem and generalize under the same problem (c.g.. by varying the grid size or PDE coefficients). For GNP achieve so, we may fine-tune a trained network with cost lower than training from scratch, or augment the GNN input with extra information (e-g., PDE coefficients). There unlikely exists an approach that remains effective beyond the problem being trained on.
Our work can be extended in many avenues. First, an immediate follow up is the preconditioning of SPD matrices. While a recent work on PDEs (Rudikov et al. 2024) shows promise on the combined use of NOs and flexible CG (Notay) 2000), is relatively fragile with respect to a large variation of the preconditioner in other problems. We speculate that split preconditioner can work more robustly and some form of autoencoder networks better serves this purpose.
Second, while GNP is trained for individual matrix in this work for simplicity, it is possible to extend it to sequence of evolving matrices, such as in sequential problems or time stepping. Con- tinual learning (Wang al. 2024), which is concerned with continuously adapting a trained neural network for slowly evolving data distributions, can amortize the initial preconditioner construction cost with more and more linear systems in sequel.
Third, the GPU memory limits the size of the matrices and the complexity of the neural networks that we can experiment with. Future work can explore multi-GPU and/or distributed training of GNNs (Kaler et al. 2022; 2023) for scaling GNP to even larger matrices. The training for large graphs typically uses neighborhood sampling to mitigate the "neighborhood explosion" problem in GNNs (Hamilton et al. 2017 Chen et al.) 2018). Some sampling approaches, such as layer-
wise sampling (Chen al. 2018). are tied to sketching the matrix product A. inside the graph convolution layer 7 with certain theoretical guarantees (Chen & Luss 2018).
Fourth, while we use the same set of hyperparameters for evaluation, nothing prevents problem- specific hyperparameter tuning when one works on an application, specially a challenging one where no purrose preconditioners work sufficiently well. We expect that this paper's results based on a naive hyperparameter setting can be quickly improved by the community. Needless to say, improving the neural network architecture can push the success of GNP much farther.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
The author would like thank eight anonymous reviewers whose comments and suggestions help significantly improve this paper.
REFERENCES
M. Arioli and I. S. Duff. Using FGMRES to obtain backward stability in mixed precision. Electron. Trans. Numer. Anal., 33:31-44, 2009.
Nathan Bell, Luke N. Olson, Jacob Schroder, and Ben Southworth. PyAMG: Algebraic multigrid solvers in python. Journal of Open Source Software, 8(87):5495, 2023.
Rishi Bommasani et al. arXiv:2108.07255 07258, 2021.
On the opportunities and risks of foundation models. Preprint
Maria Bänkestad, Jennifer Andersson, Sebastian Mair, and Jens Sjölund. Ising on the graph: Task- specific graph subsampling via the Ising model. Preprint arXiv:2402. 10206, 2024.
Raymond Hon-Fu Chan and Xiao-Qing Jin. An Introduction to Iterative Toeplitz Solvers. SIAM, 2007.
Jie Chen and Ronny Luss. Stochastic gradient descent with biased but consistent gradient estimators. Preprint arXiv: 10V7:18807118000 2018
Jie Chen and Michael L. Stein. Linear-cost covariance functions for gaussian random fields. Journal of the American Statistical Association, 118(541):177-16,, 2023.
Jie Chen, Tom L.. H. Li, and Mihai Anitescu. A parallel linear solver for multilevel Toeplitz systems with possibly several right- hand sides. Parallel Computing, 40(8):408 424, 2014.
Jie Chen, Lois Curfman McInnes, and Hong Zhang. Analysis and practical use of flexible BiCGStab. Journal of Scientific Computing, 68(2):803-825, 2016
Jie Chen, Tengfei Ma, and Cao Xiao. FastGCN: Fast learning with graph convolutional networks via importance sampling. In ICLR, 2018.
Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. Simple and deep graph convolutional networks. In ICML, 2020.
Edmond Chow and Yousef Saad. Approximate inverse preconditioners via sparse-sparse iterations. SIAM Journal on Scientific Computing, 19(3):995-1023, 1998.
Philippe Ciarlet. The Finite Element Method for Elliptic Problems. SIAM, 2002.
Timothy A. Davis and Yifan Hu. The university of Florida sparse matrix collection. ACM Trans. Math. Softw, 38(1):1-25, 2011.
S. Gerschgorin. Über die abgrenzung der eigenwerte einer matrix. Izv. Akad. Nauk. USSR Oid. Fiz.-Mat. Nauk, 6:749-754, 1931.
W. Hackbusch. A sparse matrix arithmetic based on H-matrices, part I: Introduction to H-matrices. Computing, 62(2):89-108, 1999.
W. Hackbusch and S. Börm. Data-sparse approximation by adaptive mmmmaricess Computing, 69 (1):1-35, 2002.
William L. Hamilton, Rex Ying, and Jure Leskovec. Inductive representation learning on large graphs. In NIPS, 2017.
Paul Häusner, Ozan Öktem, and Jens Sjölund. Neural incomplete factorization: learning precondi- tioners for the conjugate gradient method. Preprint arXiv:2305.16366 2023.
Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are uni- versal approximators. Newral Networks, 255):359-36,, 1989.

---

Published as a conference paper at ICLR 2025
Tim Kaler, Nickolas Stathas, Anne Ouyang, Alexandros-Statee Iliopoulos. Tao B. Schardl, Charles E. Leiserson, and Jie Chen. Accelerating training and inference of graph neural networks with fast sampling and pipelining. In MLSys, 2022.
Tim Kaler, Alexandros-Stavros Ilooouuls, Philip Murzynowski, Tao B. Schardl, Charles E. Leiser- son, and Jie Chen. Communication--fficcent graph neural networks with probabilistic neighbor- hood expansion analysis and caching. MISys, 2023.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional net- works. In ICLR, 2017.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. Journal of Machine Learning Research, 23:1-97, 2022.
Xiaoye S. Li. An overview of SuperLU: Algorithms, implementation, and user interface. ACM 325, 2005.
Trans. Mathematical Software,
Xiaoye S. Li and Meiyue Shao. A supernodal approach to incomplete LU factorization with partial pivoting. ACM Trans. Mathematical Software, 37(4), 2010.
Yichen Li, Peter Yichen Chen, Tao Du, and Wojciech Matusik. Learning preconditioner for conju- gate gradient PDE solvers. In ICML, 2023.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differ- ential equations. Preprint arXiv:2003.0348. 2020.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In ICLR, 2021.
Hong Liu, Zhiyuan Li, David Hall, Percy Liang, and Tengyu Ma. Sophia: A scalable stochastic second-order optimizer for language model pre-training. Preprint 141.22305.14342., 2023.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3:218- 229, 2021.
Thomas A. Manteuffel, John Ruge, and Ben S. Southworth. Nonsymmetric algebraic multigrid based on local approximate ideal restriction (fair). SIAM Journal on Scientific Computing, 40(6): A4105-A4130, 2018.
Thomas A. Manteuffel, Steffen Münzenmaier, John Ruge, and Ben Southworth. Nonsymmetric reduction -hase algebraic multigrid. SIAM Journal on Scientific Computing, 4165:5224- 2019.
Ronald B. Morgan. GMRES with deflated restarting. SIAM Journal on Scientific Computing, 24(1): 20-37, 2002.
M. Naumov, M. Arsaev, Castonguay, J. Cohen, J. Demouth, J. Eaton, S. Layton, N. Markovskiy, I. Reguly, N. Sakharnykh, V. Sellappan, and R. Strzodka. AmgX: A library for GPU accelerated algebraic multigrid and preconditioned iterative methods. SIAM Journal on Scientific Computing. S605)56602-666, 2015.
Yvan Notay. Flexible conjugate gradients. SIAM Journal on Scientific Computing, 22(4): 1444 1460, 2000.
M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep leamn ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
12

---

Published as a conference paper at ICLR 2025
Alexander Rudikov, Vladimir Fanaskov, Ekaterina Muravleva, Yuri M. Laevsky, and Ivan Oseledets. Neural operators meet conjugate gradients: The FCG-NO method for efficient PDE solving. In ICML, 2024.
J. W. Ruge and K. Stüben. Algebraic multigrid. In Stephen F. McCormick (ed.), Mulligrid Methods, Frontiers in Applied Mathematics, chapter 4. SIAM, 1987.
Y. Saad. IUUT: A dual threshold incomplete LU factorization. Numerical Linear Algebra with Applications, 1(4):387-400) 1994.
Youcef Saad. A flexible inner-outer preconditioned GMRES algorithm. SIAM Journal on Scientific Computing, 14(2)461-445) 1993.
Youcef Saad and Martin H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 7(3):856- 869, 1986.
Yousef Saad. Herrative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, second edition, 2003.
Vladislav Trifonov, Alexander Rudikov, Oleg Iliev, Ivan Oseledets, and Ekaterina Muravleva. Learn- ing from linear algebra: A graph neural network approach to preconditioner design for conjugate gradient solvers. Preprint urXiv:2405.15557, 2024.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In ICLR, 2018.
Liyuan Wang. Xingxing Zhang, Hang Su, and Jun Zhu. A comprehensive survey of continual learning: Theory, method and application. IEEE Transactions on Puttern Analysis and Machine Intelligence, 46(8):5362-5383, 2024.
Felix Wu, Tianyi Zhang, Amauri Holanda de Souza Jr., Christopher Fifty, Tao Yu, and Kilian Q. Weinberger. Simplifying graph convolutional networks. In ICML, 2019.
Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A comprehensive survey on graph neural networks. IEEE Transactions Neural Networks and Learning Systems, 32(1):4-24, 2021.
Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful are graph neural networks? ICLR, 2019.
Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changeheng Li, and Maosong Sun. Graph neural networks: A review of methods and applica- tions. Al Open, 1:57-81, 2020.
13

---

Published as a conference paper at ICLR 2025
A
SUPPORTING CODE
The implementation of GNP is available at whtpps//ggthhbb..om//eeeeeeeeeee GNP
B FGMRES
FGMRES is summarized in Algorithm I
Algorithm FGMRES with M being nonlinear operator
1: Let X be given. Define Hm (R 1)x+ and initialize all its entries hig to zero
2:
loop until maxi ters is reached
3: Compute rp
Axn, / = ||ro 12, and v1 = ro/B
4:
for forj= ,.... do
5:
Compute 2, = M(v,) and w = Az,
6: for i 1, ......J do
7:
Compute hij w Vi and w +- w - h;yV; end for Compute hi+ 1, w and Vj+1 w/hh+1
8:
9:
end for
Hmy 12 11: Define Zm [el...., Zm| and compute xo Zmymm where ym argminy ||āe:
12: If ||b Ax ll2 tol, exit the loop; otherwise, set Xo +- Xm
13: end loop
C PROOFS
Proof of Theorem 7
Write [|rm
S | 12 ||mm - -Fmllz It is well known that
[[ ll S K2 (X)e (m)) (
see, e.g., Saad 2003,, Proposition 6.32). the other hand, because I' = ro AZmym where y
minimizes AZmy 2a, we have Tm ro (AZm) (AZm) ) T ro Qm Q ro. Therefore,
||mm Fm 222
||Q.. Q., Qm Qm rolla & ||Qm Q Qm Qm roll2
Proof of Corollary
By taking M A - H VJ and noting that V V-1, we have
(V,Y) AM(V, Y) diag(o).......0
Then, following Theorem [rmll2 S K2 V, Y ()) We conclude the proof by noting that K2 (V, Y) (Y)..
D APPROXIMATELY EXPONENTIAL CONVERGENCE OF ||rml12
We could bound the minimax polynomial ( by using Chebyshev polynomials to obtain (ap- proximately) exponential convergence of the residual norm ||r,, for large m.
Assume that all the eigenvalues o; of Σ enclosed by ellipse E(c, d, a) which excludes the origin, where e is the center, is the focal distance, and is major semi-axis, We have
(m) Σ
min
[p(o)),
PEP
by the fact that the maximum modulus of complex analytical function is reached on the boundary of the domain. Then, we invoke Eqn (6.119) of [aad (2003) and obtain
Cm (a/dd)
(m) Σ
<
Cm (c/d)
14

---

Published as a conference paper at ICLR 2025
where Cm denotes the (complex) Chebyshev polynomial of degree m. When m is large,
Cm (maa/d) (m(c/d)
va d2 c+ v - d2
m
which gives (approximately) exponential function in Applying Corollary 2 when m is large.
conclude that
a 2 d2 c+ 62
||rmll2 2 K2 (rroll
E SCALABILITY
Based on the architecture outlined in Figure the computational cost of GNP is dominated by matrix-mand matrix multiplications, where the left matrix is either A (e.g.. in GCONV) or tall one that has rows (e.g., in the MLP and Lin layers and in GCONV). This forward analysis holds true for training and inference, because the back propagation also uses these matrix-matrix multiplications. Hence, the computational cost of GNP scales as ((nz(A) ()), the same as that of the Krylov solver.
Practical considerations, on the other hand, are more sensitive to the constants hidden inside the big- notation. For example, one needs to store many (but still a constant number of) vectors of length n, due to the batch size, the hidden dimensions of the neural network, and the automatic differentiation. Denote by è this constant number; then, the GPU faces a storage pressure nz A) + cn. Once this amount is beyond the memory capacity of the GPU, inference) by using CPU only, offloads some data to the CPU and transfers them back to the GPU on demand, undertakes more laborious engineering by distributing A and other data across multiple GPUs. For the last option (multi-GPU and/or distributed training), the literature of GNN training brings in additional techniques (such sampling and mini-batching the graph) to accelerate the computation (Kaler et al. 2022 2023).
either performs training (and even
F ADDITIONAL. EXPERIMENT RESULTS
F.1 ADDITIONAL METRICS FOR COMPARING PRECONDITIONERS
In addition to the Iter-AUC and Time-AUC metrics, one may be interested in using the final relative residual norm compare the preconditioners. To this end, summarize the results in Figure by using this residual metric, as an extension of Figure 2 Note that for each problem, we solve the linear system twice, using rtol and max ters the stopping criteria; and the other time using rt o 1 and t timeout. Hence, there are two charts in Figure 7 depending on how the itcrations terminated.
by residual terminated by maxiters
by residual
terminated by timeout
preconditioner
37.7%
16.0%
37.3%
EU
CNP
0.3%
1.2%
Jacobi
14.5%
2.2%
29.2%
13.8%
GMRES precond
No precand
Figure 7: Percentage problems on which each preconditioner performs the best.
We see that across these two termination scenarios, observations similar. Compared to the use of the Iter-AUC and Time-AUC metrics, GNP is the best for a similar portion of problems under the residual metric. The portion of problems that ILU and GMRES (as preconditioner) perform the best increases under this metric, while the portion that AMG performs the best decreases. Jacobi and preconditioner remain less competitive.
15

---

Published as a conference paper at ICLR 2025
F.2 PRECONDITIONER CONSTRUCTION TIME AND SOLUTION TIME
In Figure S we plot the preconditioner construction time and solution time for each matrix and each preconditioner. This figure is counterpart of Figure 3 by removing the log scale in the y axis. One sees that the GNP construction time nearly proportional to the matrix size.
Thmee (seconds)
300
7500
0.75
100
0.50
10
100
50
2500
0.25
0.00
10
10
10
10
matrix
Construction
Solution
Figure 8: Preconditioner construction time and solution time (using maxiters to stop). The axis the linear scale.
F.3 BREAKDOWN OF BEST PRECONDITIONERS WITH RESPECT TO MATRIX DISTRIBUTIONS
In Figure 9 we plot the proportion of problems on which each preconditioner performs the best, for different matrix sizes, condition numbers, and application areas. This figure is a counterpart of Figure 4 The main observation that the proportion of problems on which GNP performs the best is relatively flat across matrix sizes. Other observations follow closely from those of Figure 4 namely, GNP particularly useful for ill-conditioned problems and some application areas (chemi- cal process simulation problems, economic problems, and eigenvalue/model reduction problems).
Beet precerdiiioee
08
wasls Witsrution
proportion
AMG
0.6
04
0.4
Jacobi
GMHESS as procond
pracond
0.0
10
matr size 1C:04
leg10icondition number)
Proportion of problems on which each preconditioner performs the best, for different matrix sizes, condition numbers, and application areas. Only the application areas with the top number of problems are shown. The last bar in the middle plot is for condition number 2 1026 Note that not every matrix has a known condition number.
Figure 9:
F.4 COMPARISON BETWEEN PYAMG AND AMGX
Besides PyAMG, another notable implementation AMG (which supports using as precondi- tioner and supports GPU) is AmgX / thub. com/NVTDDA/ANT AMMGX Here, we com- pare the two implementations in terms of robustness and speed.
Unlike PyAMG, which is blackbox, AmgX leaves the choice of the AMG method and the parameters to the user. We used the driver examples examples/amgx.ccage /aagxxccapp....... a and the configuration sre/conficn FGGRRS AGGREGATION json, which appears be the most reasonable for the diverse matrices and application areas at hand. This setting uses aggregation- based AMG rather than classical AMG, which is also consistent with the preference expressed PyAMG.
Here are some findings.
AmgX is less robust than PyAMG. Table 3 shows that AmgX causes more failures during the struction as well as the solution phase. In fact, there were three cases where AmgX hanged, which hampered automated benchmarking.
Despite a lack of robustness, the construction of AmgX is much faster than that of PyAMG. For AmgX, the distribution of the preconditioner construction time has minimum 0.006s, median
16

---

Published as a conference paper at ICLR 2025
Table 3: Failures two implementations of AMG.
PyAMG
AmgX
Construction failure
62 (7.15%)
105 (12.11%)
Solution failure
(0.58%)
36 (4.15%)
0.017s, and maximum 4.895s. In contrast, the distribution for PyAMG has minimum 0.005s, median 0.286, and maximum 8505. Is. Nevertheless, there are still 16% of the cases where the PyAMG preconditioner is faster to construct than the AmgX preconditioner.
We skip the comparison of the solution time. This is because in AmgX, the Krylov solver and the preconditioner are bundled inside a C complementation. A fair comparison requires isolating preconditioner from the solver, because our benchmarking framework implements the solver in Python. However, the isolation is beyond scope due to the complexity of the implementation, such as the use of data structures and cuda handling. Nevertheless, we expect that the solution time of AmgX is faster than that of PyAMG, based on the observations made for preconditioner construction.
F.5 CHANGING THE DEFAULT PYAMG SOLVER
The blackbox solver of PyAMG, pyamg, pyamg.blackkal..........k blackbox solver, uses the classical-style smoothed aggregation method pyang aggregation smoothed aggregation solver (SA). For nonsymmetric matrices, however, the AIR method works better (Manteuffel et al.) 2019 2018). In this experiment, we modify the blackbox solver by calling pyamg. clwss ical. Lr_solver for nonsymmetric matrices.
We first tried using the default options of AIR, but encountered "zero diagonal encountered cobi" errors. These errors were caused by the use of jc_jacobi as the post-smoother. Then, we replaced the smoothers with the default ones used by SA (gause_seidelng) solver hung the problem FIDAP exx40. Before this, 105 problems were solved, among which 83 were nonsymmetric. The results of these 83 systems are summarized in Table suggesting that AIR is indeed better but not always.
In this case, the
Table 4: Counts of problems when comparing AIR with SA for nonsymmetric matrices.
AIR is is better than SA
AIR is worse
The two are equal
By iteration count By final residual
14
50
29
Overall, we can conclude that AIR has not exhibited its full potential for nonsymmetric matrices. We speculate that robust implementation of AIR challenging and that might be the reason why the authors did not use it for the blackbox solver in the first place.
F.6 CONVERGENCE HISTORIES
In Figure 10 we show the convergence histories of the systems appearing in Figure 5 with re- spect to both iteration count and time. We also plot the training curves. The short solution time of GNP is notable, especially when compared with the solution time of GMRES as a precondi- tioner. The training loss generally stays on the level of 1e-2 1e-3. For some problems (e.g., VanVel he.eellee_std___p0
training does not seem to progress well, but the training loss is already satisfactory point in the first place, and hence the preconditioner can be rather effective. We speculate that this behavior is caused by the favorable architectural inductive bias in the GNN.
To provide further examples, areas identified by Figure 4 to favor GNP: chemical process simulation problems (Figures II 13), economic problems (Figures 16), and eigenvaluc/model reduction problems (Figures) 17 18). In many occasions, when GNP performs the best, it outperforms the second best substantially. It can even decrease the relative residual norm under rtol while other preconditioners barely solve the system.
show the convergence histories of all systems in the application
17

---

Published as a conference paper at ICLR 2025
GHS ji_jnnd/e/b/ddddd..mm 2D/3D problem, cond 5.92e+02
GNP
log10(relative residual)
reelltieee residual)
-0.75
AMG
Jacobi
log10(training
-1.00
GMRES as precond
No precond
-2
-6
-1.75
100
1500
2000
iterations
time (seconds)
epochs
VanVelzen/stdl Jac3, chemical process simulation problem, cond 1.02e+27
GNP
log10 residual)
log10 residual)
logl0(training loss)
-2.825
AMG
Jacobi
-2.850
GMRES as precond
-2.875
No precond
-2.900
-2.925
-8
-8
-2.550
100
500
1000
1500
2000
iterations
time (seconds)
epochs
Simon/venkat25, subsequent computational fluid dynamics problem, cond = N/A
GNP
log10(relative residual)
(reeatvee residual)
log 10trainiig loss)
AMG
-2.6
Jacobi
-1-
-1
GMRES as precond
-2.8
precond
-2
-3.0
25
50
75
100
2
500
1000
1500
2000
iterations
time (seconds)
epochs
Hollinger/mark3jac0000, economic problem, cond 6.79e+13
GNP
log orrellaeeevve residual)
log1Ofrelative residual)
19g1ttraining loss)
AMG
-2.50
Jacobi
-1-
GMRES as precond
-2.55
No precond
-2-
-2.60
-3-
-2.65
25
500
1000
2000
iterations
time (seconds)
epochs
Figure 10: Convergence the linear system solutions and training history of the preconditioners,
for the matrices shown in Figure 5
18

---

Published as a conference paper at ICLR 2025
wB/test1505 cand L.59e+ 1/
HB/vest2021
Baibbrr2000
Grundib22 ss cond /.bb++11
cand 1.50e
cond 1.
log logelireee residual)
-1
-5
-2
-2-
-3-
-3-
-10-
-4
-4
-4
-13
0
50
100
0
50
100
0
50
100
50
100
Grund/h dyn
Grund/bayer01
Brund/bayer02
Grund[tayer03 cond 1.59c+19
cond 2.53e+ 14
5,07e+10
log10/eelative resicual)
-2
-4
-4-
-6-
-6-
-5
-8
-3
30
100
30
100
50
100
50
100
Grund/bayer04
Crund/bayer05 1.96e-30
Grund/beyer06
Grund/bayer07
1.60e+28
3.79e429
log1ofrelative reciduall
0.0
0.0
-2.
-5.0
-4
-5.0
-3
-7.5
-4
-10.0
-3
-10.0
50
100
50
100
50
100
50
100
Grund/aayeri88
Grund/bayer09
Grund/bayer30
Mallya/hhr0
5,57e+27
6.33e+05
logi0/relstion res dual)
-2-
-2
-4-
-2
-6
6-
-8
-8
50
100
50
100
50
100
50
100
Mallya/hrr2
Mallya/hhr4
Mallys/(hrO4c
Mallya/hrr07
4.12e-118
3.80c
-4.40c
log10/rel ative recidual)
0.000
0.0
0.0
-2
-0.025
-0.2
-4
-0.030
-0.4
-0.075
-0.2
-0.100
-0.3
0..0
50
100
50
100
50
100
50
100
iterations
AMG
Jacobi
precond
No precond
Figure 11: Convergence of the linear system solutions for chemical process simulation problems (1/3).
19

---

Published as a conference paper at ICLR 2025
Mallys//hr07c
Mallya/hr10
Mallya/thr10c
Mallya/th/11
cond 2.39e+13
cond =1.33e+18
cond 7.39e+13
cond 6.42e+19
log10 lative resiuall)
0.000
0.00
0.00
-0.025
-0.25
-3.5
-0.01
0..00
-0.02
-0.75
-0.075
-1.00
-1.5
0.03
50
100
0
50
100
50
100
50
100
Mal
Mallya/hr14 cond 8.57e++22
Mallya/hhr24:
Mallya/thr17
cond 3.96e+14
cond 4.03e+14
cond 4.22e+23
log10/relatiee resicual)
0..00
3.0
0.00
0.0
-0.05
-0.02
-3.5
-0.10
-0.04
-10-
-1.0
-0.15
0
50
100
0
50
100
50
100
50
100
Mallya/thr17c
Mallya/hrr4
Mallya/Ihr334c
Mallya/th
3,85p
4.
cond 1.80e 116
logtt/relattv residuall
0.0
3.
000
0.0
-0.1
0..02
-0.5
-0.2
-1.0
-0.04
50
100
50
100
50
100
50
100
Mallye/thTTcc
Zimee/ettr
Ztneyhydri.
Zitneyiradr11
3.377+12
-2335e+14
5.93e+10
loglorreiative res dusidaal)
0.0
0
0.0
-0.1
-2-
-3.0
-7.5
-0.3
-3
100.0
50
100
50
100
100
50
100
Zitnev/diist1
Zihney/cist2
Zineerddistaa
Venveeren/stt1 Jaz7
1.74p
1.3774
-1.2e
log1trrai lativa residual)
0.0
-2
0.5
-5
50
100
50
100
100
50
100
iterations
precond
precond
Figure 12: Convergence of the linear system solutions for chemical process simulation problems (2/3).
20

---

Published as a conference paper at ICLR 2025
VarVVelaen/std__eek Jac2
VanVelzen/stol__ert Jac3_ db cond 4.15e+22
Ven//eiee/ttd1__B 1acc cond 1.02+277
VanVelzenizd Ja:2 db cond 5.95e+18
cond 1.88e+26
log1 Direraativve residual)
0.00
-0.25
-0.50
-6
-6
-1 11.00-
100
103
50
VanVelzervzd Jac2
VarvelzennVZ Jac3 db
VanVelzenvZZ JaCS cond 2.45e+24
VanVelzen/2d
cond 1.95e+224
cond 5.59e+18
cond 2.72e+18
log DgllDieelatve residua))
0.0
-0.2
-2
-02
00.4
-4
-0.5
-4
-0.4
30
100
0
30
103
50
100
30
100
VanVelzenvZd Jac6
cond 1.53e+26
log10irelsten lativee res dual)
-2-
-4-
100
Iterations
Jocoii
precond
No precond
Figure 13: Convergence of the linear system solutions for chemical process simulation problems (3/3).
21

---

Published as a conference paper at ICLR 2025
1B//nahiniaa cond 2,08e+13 13
HB/arani678 cond = 9.50te+U
HB/psmigr 1 cond 842+39
HD/psmigr 2 cond L./2e+04
log10irelative residual)
-2
-2
-2
-6-
-6
-8%
-8
-8
50
100
50
100
50
100
100
HB/psmigr gr 3
Grund/poli cond 3.12a
Grund/poli large cond NAA
Hallingerg7/acc100
cond 1.33c.02
cond 1.03e 1.
loglovelltive resiuaal)
3.0
0.0
-2.5
-25
-5.0
5500
-0.4
-10-
-7.5
-7.5
-15
-10.0
0
10
20
30
0
25
30
75
0
50
100
50
100
7jlingerr/7/ac0000
Hollingerig7 ac020sc
Hallingerlg]JacC40 cund 1.45++17
Lond 3.50e+13
coral 5,26e=16
cond 4.45e+12
log10irelative res duall
0.00
5.0
0.0
-0.25
-0.50
-0.5
-0.5
-1.0
-1.0
-3-
0
50
100
0
50
100
50
100
50
100
Holligerig7/ac050s
Hollinger/g7jac000 cond 5.ee+117
Hollinge-7//aa000006 cond 6.74e+12
4.34e+12
cond 9.35e+112
1eg residuall
0.
0.0
-0.5-
-0.5
11
-0.5
-2-
-1.0
-1.0
-1.0
-1.5%
0
50
100
0
50
100
50
100
50
100
Hollinoera77:a0000
Hollingerig/aac00ss
Hollinger/o7aac1000
Hollinge cond 8.57e+
2.41e+18
6.80e+12
6.87e
teglorelative residual)
0
0
0.0
-2-
-05
-05
-1.0
-4-
-1.5
-1.5
0
50
100
0
30
100
50
100
50
100
iterations
ILUI
AMG
jocobi
GMRES
precond
No precond
Figure 14: Convergence of the linear system solutions for economic problems (1/3).
22

---

Published as a conference paper at ICLR 2025
Holliggero7|ec120 cond 2.69e+19
Holl ngerlg7|ac1205cc
Hellingert7/jac445 cond N/A
Hbllingerig/aac14DSS
cond = 1.01e+13
cond = N/A
log10/relative residual)
-0.5
-2
-0.5
-2
-3-
-1.0-
-1.0
-4
50
100
100
100
100
Hollingerig/jec100 cond NJA
Holl ngerigg/jc180000 cond - M/A
Hallingergg/jac1ek cond N/A
H
cond N/A
log20/eelatve resiuall)
0.0
0.0
-0.5
-2
-1.0
-1.0
-6
30
100
50
100
0
50
100
50
100
Hollingertg/iac200
Holl nger/g7/ac200ssc
Hollingerian9jjac220
Hbllinger/an39jaac200s cond 2.08e-00
cund
N/A
cond 3.19e+13
log10irelatic res duall
0.0
-0.
-2-
-0.5
-1.0
4
-1.5
-4-
-2.0
50
100
50
100
50
100
50
100
Hollinger/jn99jjc040
Hellings//jnn@jaa040o:: cond 4,29e+000 06
Hollinger(an@@joc066 cond 1.33e+14
Hbllingerjjn9ijxx000c cond 6.77e-06
cond 7.47e+13
le residuall
0.0
-1
-1
0.5
-2-
-1.0
-3
50
100
50
100
0
50
100
50
100
Hollnngrrjan99ac0B0
Hollinge:/ann@@aa00000 cand 8.23e D6
#ollingertan99jaa1000
Hollinger/jan?]jac
cond 2.33.c
3.79c 114
cand 9.82c 06
logt0irelative residual)
0.0
0.0
-1-
-0.5
-2-
-1.0
-1.0
-1.5
50
100
100
50
100
100
Iterations
AMG
jocobi
GMRES
precond
No precond
Figure 15: Convergence of the linear system solutions for economic problems (2/3).
23

---

Published as a conference paper at ICLR 2025
Holl nger)an99jac1100 cand NVA
Helllnggrrarn9gjaa1200 cond I/A
#bllinger/marr3(@c220 cond = 3.E9e+13
Hollingerimark/jac0a cond = 8.40e+13
log10/nelative residual)
-0.25
1
-0.50
-2-
-1.00
-3
-3
50
100
50
100
50
100
50
100
#bllingger/mar33joc000 cond 6.79c4113
Hollingermmakkjppclke cond 5.51a
Hbllingermmrkkjac050 cond 1.13pi
Hollingermark
cond - 5,50e 13
leg20(eelaive residual)
-1-
-1
-2
-2
-3
50
100
0
50
100
0
50
100
0
50
100
#bllinger/mark3jjcc000
Hollinger/mark3jact00ss
#bllingerimark3/ac1100 Lond NVA
Hollinger/markjjac100a
cond 2.c8e+14
cond 5.630+13
cund N/A
log10irelative dusi)
0.0
-0.5
-1
-0.5
-1.0
-1.0
-15
-15
-2.0
-2.0
-2
0
50
100
0
50
100
0
50
100
0
50
100
#bllinger/mark3jcc220
Hbllingge/mmrt3jaa1220e cond N/A
Hbllinger/markjje2110
Hollinger/mar3jjc114.00
zond N/A
cond N/A
cord N/A
ve residuall
0.00
0.0
-.2
-0.50
0.4
-0.55
-2
-1.00
0
50
100
0
50
100
0
50
100
0
50
100
Marini/eurqsa
Grundrpcli3
Grundrodli4
2.81e120
cond S.25c1100 02
cond 7.50c1
logl0relatiiv residual)
0.0
-2
-1-
-5.0
-2-
-7.5
-10.0
0
50 Iterattonn
100
0
30
100
50
100
AMG
jocobi
GMRES precond
No precond
Figure 16: Convergence of the linear system solutions for economic problems (3/3).
24

---

Published as a conference paper at ICLR 2025
Rommes/www vrof 6405
Rommes//ingo, afoossc Itaip cord 0.54e+445
Rommes/mimoars sysstmm /.1be
Rommes//mimo2Bii2 system
b3c+444
cand 4.84444
log residual)
-5
-5
-10-
-10
-10
-15
-15
2
0
Rommes/m imo46x46 system 77.39e+29
Rommesjjuba00K
Rommes/bauru5727
Rommesizzrrs,, ropss 13k cand 2.34c+55
logldrrelative resicual)
0.0
0.0
-5
-2.
-2.5
-10
-5.0
-10
-7.5
-15
-15
10.0
-10.0
50
100
30
100
0
Rommesidescriptor xingo5u 6.77e43
Rommesinopss 11k
Rommessxingg0112
Rommes/b ps98 606 5.8le+26
050++65
1.777+33
logiolnelative reciduall
0.0
0.0
-5
-5.0
-10
5.00
-7.5
-15
77.5
-15-
10.00
10.00
-20
50
100
2
50
100
Rorrmes/bins96. 1142
Rommes/bipss8 1450
Rommes/bipsd7, 1199
Rormees/biss77_ 1998
9.95e+32
2.05++30
3.10e+26
=2.99=229
logt0irelative desidus))
-5-
-5
-10
-10
-5-
-15
-15
-10
0
2
50
100
Rorrmes/bias07 2475
Rommes/bips07,
Rammes/M20PI
Rommes/M40PI
1.66ci.36
log10/rel stive residual)
0.0
0.0
-2.5
0.5
-5.0
-1.0
-7.5
-10
-1.5
-1.00
-2.0
50 iterations
100
0
50
100
0
50
100
50
100
precond
precond
Figure 17: Convergence of the linear system solutions for eigenvalue/model reduction problems (1/2).
25

---

Published as a conference paper at ICLR 2025
Rommas/M20PP Ln cond 6.00e+177
Rommes/M10P1 cond 3.93e+18
Rommes/MM8PL
cond 1.38e+10
cond 3.53e-17
Drrrraaaevvee residual)
0.0
0.0
0.0
0.0
0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1.5
-1.5
100
100
100
Rommes/520PI
Rommes/544P)
Rommes/580P)
Rommes/520FL
cond 1.41e+110
cond 4.80e
cond 7.38e
cond 6.b0e-17
lgg residaa))
0.0
0.0
00.5
-0.5
-0.5
-1.0
-1.5
-1.0
-2.0
30
100
30
100
0
30
100
0
30
100
Rommes/SAoPL
Rommes/SSB0PIn cond 3,63e+17
cond 3.35e+18
log1 Direlativv residual)
0.0
0.0
-0.3
-0.5
-1.0
-1.0
-1.5
100
100
iterations
AMG
Jocobi
precond
No precond
Figure 18: Convergence of the linear system solutions for cigenvalue/model reduction problems (2/2).
26


Reviews:
[{"review_id": "ps8G0xTJHx", "reviewer": "Reviewer_VMvP", "summary": "Authors propose use Graph Neural Network (GNN) as a general-purpose preconditioner. They offer a convergence analysis for flexible GMRES. A new effective training data generation is proposed to training GNN. They develop a scale-equivariant GNN as the preconditioner. A novel evaluation protocol proposed by authors is being for general-purpose preconditioners.", "questions": "See weaknesses.", "limitations": "", "rating": 6, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "* Authors pay their attention to the problem which is linked to using neural networks (in particular, GNN) as the preconditioner. In this sense, the preconditioner is a nonlinear operator and, using a standard preconditioned solver with this preconditioner is not correct. In this case, the convergence theory is broken. They concentrate on flexible GMRES and present an original convergence analysis. \n* Authors propose original scheme of data generation. Suggested sampling is linked with eigen-subspace of the $\\mathbf{A}$.\n* Of course, the idea of using GNN as the preconditioner is not novel. Authors offer a fresh approach to normalize $\\mathbf{A}$. This prevents the potential division-by-zero issue that can arise in the standard normalization of GCN. The major innovation of the architecture is its scale-equivariance, i.e. the input space of the neural network is restricted.\n* Proposed architecture is checked on the SuiteSparse matrix collection https://sparse.tamu.edu (square, real-valued 867 matrices (not spd-matrices) from 50 application areas, whose count of non-zero elements is less than 2M. To compare with classical methods such as ILU, AMG, and GMRES author propose two novel metrics.", "weaknesses": "* Authors use $\\ell_1$ residual norm as the training loss. Is it possible to use $\\ell_2$ norm or others as the training loss?\n* Authors compare their architecture with classical methods(GNP, ILU, AMG, Jacobi, and GMREs). It will be good to compare their approach and metrics with the existing general-purpose preconditioners which are used GNN as the preconditioner. See for example, https://proceedings.mlr.press/v202/li23e/li23e.pdf, https://arxiv.org/pdf/2405.15557.\n* Are there exist some constraints to use this approach for SPD-matrices?\n* Why is not the full GitHub repository made available for review?\n* Why do the authors only assume the ground truth solution $\\mathbf{x} = 1$?", "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": "", "extra_scores": {}}, {"review_id": "9zHZ4FDbDf", "reviewer": "Reviewer_7HSM", "summary": "The paper proposes a novel approach using Graph Neural Networks (GNNs) as general-purpose preconditioners for solving large, sparse linear systems. The authors present an empirical evaluation of their method, showing that it outperforms traditional preconditioning techniques like ILU and AMG in terms of construction time and execution efficiency, especially for ill-conditioned matrices. The method is tested on a diverse set of over 800 matrices, suggesting significant advantages in robustness and performance compared to existing methods.", "questions": "1. How does the proposed preconditioner handle matrices outside the tested conditions, particularly in real-world applications?\n\n2. Could you elaborate on the training data generation process? How might biases in the training set affect the performance of the GNN?\n\n3. What are the implications of your findings for future research in preconditioning techniques, especially concerning adapting the GNN for a wider array of problems?", "limitations": "", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- **Innovative Approach:** The use of GNNs as preconditioners is a fresh perspective that leverages recent advancements in machine learning.\n\n- **Empirical Validation:** The extensive evaluation on a wide range of matrices provides strong evidence of the proposed method's effectiveness across various problem domains.\n\n- **Robustness:** The proposed method demonstrates high robustness with a low failure rate compared to traditional preconditioners, particularly in handling ill-conditioned problems.\n\n- **Theoretical Contributions:** The paper includes a convergence analysis for the flexible GMRES method, which adds valuable theoretical insights to the practical implementation.", "weaknesses": "- **Scalability Concerns:** While the paper mentions robustness, it lacks detailed discussions on scalability, particularly regarding the performance of GNNs with significantly larger matrices.\n\n- **Dependence on Training Data:** The effectiveness of the GNN-based preconditioner seems to heavily rely on the quality and diversity of the training data. The paper does not explore the implications of this dependence thoroughly.\n\n- **Limited Theoretical Foundation:** Although the empirical results are compelling, the theoretical grounding could be strengthened, especially concerning the convergence properties of the proposed methods under varied conditions.", "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": "", "extra_scores": {}}, {"review_id": "tpUWOOkVqF", "reviewer": "Reviewer_HJDk", "summary": "This paper proposes a GNN-based approach to preconditioning Krylov solvers, offering a novel alternative to traditional algebraic preconditioners. Unlike conventional methods, the GNN-based preconditioner (GNP) approximates the matrix inverse without relying on additional information, making it adaptable and robust across a wide range of matrices. GNP demonstrates predictable training costs, good robustness, and, in some cases, faster convergence than traditional preconditioners like ILU and AMG, which are often hampered by unpredictable construction times and structural limitations.\nThe contributions:\n1. Convergence analysis for FGMRES, a widely used but theoretically underexplored method.  \n2. New approach to training the neural preconditioner, focusing on effective training data generation by sampling from the bottom eigensubspace of  A to enhance training performance. \n3. Scale-equivariant GNN to serve as the preconditioner, addressing the challenge of varying input scales in the training data by enforcing an inductive bias that maintains scale-equivariance. \n4. Novel evaluation protocol to test the preconditioner broadly, evaluating across over 800 matrices from 50 diverse application areas in the SuiteSparse collection.", "questions": "1. Since neither Iter-AUC nor Time-AUC specifically prioritizes the final residual accuracy, have you considered including a metric that directly measures the final residual? For applications where achieving a specific residual threshold is crucial, this could provide a more balanced evaluation. If not, what are your thoughts on this?\n2. Iter-AUC and Time-AUC capture different aspects of convergence (iteration efficiency and time efficiency, respectively). In cases where these metrics might conflict (e.g., one method scores well on Iter-AUC but poorly on Time-AUC), how would you recommend interpreting the results?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1. Introduces a GNN-based preconditioner (GNP) as a novel, general-purpose alternative to traditional methods, approximating the matrix inverse without additional problem-specific information.\n2. Offers stable and consistent training times, unlike traditional preconditioners with unpredictable construction times due to matrix irregularities.\n3. Contributes a convergence analysis for FGMRES.\n4. Tests GNP on over 800 matrices from 50 application areas.", "weaknesses": "Sometimes both Iter-AUC and Time-AUC do not give good \"weight\" for the final residual accuracy. A method that reaches a lower final residual can be underrepresented if it converges more gradually, while a method with fast early convergence but a higher final residual might appear more favorable.", "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": "", "extra_scores": {}}, {"review_id": "0yF0aMFiPq", "reviewer": "Reviewer_4sht", "summary": "The authors focus on the task of preconditioning linear systems which remains one of the most important problems in scientific computing and applied maths. The authors use a machine learning approach and propose using graph neural networks as a general-purpose preconditioner.", "questions": "- For the training data generation at the beginning of 2.2, is this really meaningful? You cannot sample the whole of $\\mathbb{R}^n$ to run the problem forwards as also given possible ill-conditioning of $A$ approximating the inverse is difficult. Is this the reason the authors then switch to Arnoldi? This is the basis for many Krylov solvers and it is not clear what the advantage of the machine learning approach his if we have to do Arnoldi anyways?\n- How about the sharpness of the Gershgorin estimates? These are not necessarily sharp.\n- Do the authors really observe convergence of GMRES in 10 iterations to $10^{-6}$?\n- AMG works better for symmetric problems, for nonsymmetric matrices PyAMG comes with better approaches for nonsymmetric problems such as the AIR approach.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Preconditioning is a relevant problem and discovering new preconditioners for novel applications is resource-demanding\n- The paper tackles a large number of test matrices and compares with SOTA black box preconditioners", "weaknesses": "- The authors state their result as probably novel. Have they not found anything similar in the literature? I would like them to discuss the difference to the work USING FGMRES TO OBTAIN BACKWARD STABILITY IN MIXED PRECISION by Arioli and Duff, which they have done in the revised version.\n- I am not a fan of presenting the results as best performer but it could be that the other methods are really close and this method only outperforms slightly and is terrible for the other problems. I would like to see actual  numbers or convergence plots.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)