# Towards Neural Sparse Linear Solvers

Luca Grementieri  Paolo Galeone

###### Abstract

Large sparse symmetric linear systems appear in several branches of science
and engineering thanks to the widespread use of the finite element method
(FEM). The fastest sparse linear solvers available implement hybrid iterative
methods. These methods are based on heuristic algorithms to permute rows and
columns or find a preconditioner matrix. In addition, they are inherently
sequential, making them unable to leverage the GPU processing power entirely.
We propose neural sparse linear solvers, a deep learning framework to learn
approximate solvers for sparse symmetric linear systems. Our method relies on
representing a sparse symmetric linear system as an undirected weighted graph.
Such graph representation is inherently permutation-equivariant and scale-
invariant, and it can become the input to a graph neural network trained to
regress the solution. We test neural sparse linear solvers on static linear
analysis problems from structural engineering. Our method is less accurate
than classic algorithms, but it is hardware-independent, fast on GPUs, and
applicable to generic sparse symmetric systems without any additional
hypothesis. Although many limitations remain, this study shows a general
approach to tackle problems involving sparse symmetric matrices using graph
neural networks.

sparse linear solvers, graph neural network, numerical analysis,

##  1 Introduction

The resolution of linear systems is a classical problem at the core of the
numerical analysis. Linear systems are so pervasive because they are a
fundamental component of the numerical resolution of differential equations.
For this reason, they appear in all branches of science: from engineering to
biology. The finite element method (FEM) is a popular and accurate method to
solve differential equations numerically (Logan, 2016). In this context, the
matrix associated with the linear system has the property to be very sparse
and symmetric (and often positive-definite). The real-time simulation of FEM
models requires fast resolution algorithms that can be inaccurate up to a
certain degree. Classic methods to solve sparse linear systems have been
conceived before the introduction of GPGPU computing when the processing power
available was very limited, so they are not suited to use a GPU effectively.
The methods to solve sparse linear systems are classified into direct and
iterative methods.

Direct methods, also known as decomposition methods, rely on a factorization
of the coefficient matrix in a pair of matrices that have simplifying
properties, like triangular or orthogonal matrices. The most established
direct methods are LU factorization (equivalent to Cholesky decomposition for
definite-positive matrices) and QR decomposition. Decomposition methods are
very accurate, but their parallelization on modern hardware (GPUs) is
complicated, and it is impossible to trade off accuracy for computational
time. Moreover, the factors obtained in the decomposition process can be much
less sparse than the starting matrix, causing an increase in memory
consumption. Every null element of the initial matrix that becomes a non-null
value during the factorization is called a fill-in. When the number of fill-
ins is large, the execution time increases considerably because the mutation
of the sparsity pattern of a sparse matrix is a slow operation. A permutation
of rows and columns can mitigate the appearance of fill-ins by concentrating
the non-null elements of the sparse matrix around the diagonal. The research
of the optimal permutation that minimizes the number of fill-ins is an NP-hard
problem (Yannakakis, 1981), so heuristic fill-reducing algorithms are employed
in practice.

Iterative methods repeatedly refine an approximate solution starting from an
initial guess, often the null vector. Iterative methods have the advantage
that a less stringent stopping criterion permits the reduction of the
execution time at the expense of a less accurate solution. The coefficient
matrix is never changed or inverted, but it is used only to perform matrix-
vector multiplications. For this reason, iterations are fast on GPU but, since
these methods are inherently sequential, the GPU usage is low, and the degree
of parallelization is limited. When the condition number of the coefficient
matrix is high, iterative methods can incur a slow convergence increasing the
total computation time. Preconditioning methods are harnessed to reduce the
condition number, but they can be slow on GPUs because they often rely on
decomposition methods. The most adopted iterative methods are Conjugate
Gradient (CG) or its variant Preconditioned Conjugate Gradient (PCG). Such
algorithms are only applicable to positive-definite matrices. When the
coefficient matrix is indefinite, it is possible to employ BiCGSTAB (Van der
Vorst, 1992) or GMRES (Saad & Schultz, 1986), but they are often slower than
CG on systems of similar size.

We are interested in fast approximate solvers that can be used in real-time
non-critical applications to predict a coarse solution of a linear system. We
present a deep learning framework to learn approximate sparse linear solvers
tailored for a specific use case. Our approximate solvers hinge on the recent
advances in graph neural networks (GNNs) and geometric deep learning
(Bronstein et al., 2017). To leverage GNNs, we present a novel way to
represent a sparse symmetric linear system as a weighted undirected graph.
Using this input representation, we can cast the resolution of a sparse linear
system as a node regression task. We tackle this regression task using a graph
convolution network (GCN) coupled with some techniques to exploit the scale
invariance property of the problem.

We test our neural sparse linear solvers (NSLSs) on a dataset of synthetic and
realistic symmetric positive-definite linear systems coming from the field of
structural engineering. Our tests show that our approximate solver is less
accurate than existing methods and, for this reason, slower if compared with
other iterative methods (being the relative error equal). Nevertheless, an
NSLS has several advantages: it does not suffer from convergence issues caused
by ill-conditioned matrices; its implementation is hardware-independent; the
training on positive-definite, indefinite, or general sparse symmetric systems
does not require any change.

Lastly, we consider our proposed deep learning approach an initial significant
example of how deep learning can help solve numerical analysis problems on
sparse matrices. Our technique is part of the broader line of research toward
computational-aided mathematics (Lample & Charton, 2019; Davies et al., 2021).

To summarize, our contributions are the following:

  * •

we introduce a novel representation of a sparse symmetric linear system as a
weighted undirected graph;

  * •

we frame the resolution of a sparse linear system as a node regression task;

  * •

we build a GNN-based deep learning model to solve approximately linear
systems.

##  2 Background

###  2.1 Problem Definition

Let $Sym_{n}$ be the set of real symmetric $n\times n$ matrices. A symmetric
linear system is defined by a symmetric matrix of coefficients $A\in Sym_{n}$
and the vector of constant terms $\mathbf{b}\in\mathbb{R}^{n}$. The solution
of the linear system is the vector $\mathbf{x}\in\mathbb{R}^{n}$ such that

$A\mathbf{x}=\mathbf{b}.$ (1)

A system is said sparse if $A$ is a sparse matrix, i.e. the number of non-zero
entries scales linearly with the number of rows and columns.

If $\det(A)\neq 0$, the linear system (1) has a unique solution, so it exists
a solution function $s$ that maps the linear system to its solution

$s(A,\mathbf{b})=A^{-1}\mathbf{b}=\mathbf{x}.$ (2)

$s$ is linear with respect to $\mathbf{b}$, but it is non-linear with respect
to $A$ because matrix inversion is a non-linear operation.

The computation of $s$ is slow when $n$ is large, so we look for an algorithm
to approximate the function $s$.

###  2.2 Preconditioning

One possibility is to approximate the inversion of $A$ separating the linear
part of $s$, the multiplication with $\mathbf{b}$ from the non-linear
inversion of $A$. The line of research that proposes new preconditioning
methods follows this approach implicitly.

A preconditioner $M$ of a matrix $A$ is a matrix such that $M^{-1}A$ has a
smaller condition number than $A$ (and $M^{-1}$ is easy to compute). When a
preconditioner $M$ is available, the solution $\mathbf{x}$ is computed solving
the equivalent system $M^{-1}A\mathbf{x}=M^{-1}\mathbf{b}$. For a symmetric
matrix $A$ (a special case of a normal matrix), the condition number is

$\kappa(A)=||A^{-1}||\cdot||A||=\frac{|\lambda_{\text{max}}(A)|}{|\lambda_{\text{min}}({A})|}\geq
1,$

where $\lambda_{\text{max}}(A)$ and $\lambda_{\text{min}}(A)$ are maximal and
minimal (by modulus) eigenvalues of $A$ respectively. It follows that the
minimum condition number is attained on multiples of the identity matrix.

Preconditioning methods try to approximate the inverse of $A$ under some
restrictive conditions to minimize the condition number. For example, usually,
the sparsity pattern of $M$ or $M^{-1}$ coincides with the sparsity pattern of
$A$. With such limitation, the inverse $A^{-1}$ can be approximated but not
recovered because the inverse matrix of a sparse matrix is non-sparse in the
general case. The methods that learn the preconditioning method based on a
dataset of matrices often impose the same sparsity constraint (Götz & Anzt,
2018).

##  3 Graph Representation of a Sparse System

In this work, we propose a method to approximate the function $s$ directly
using a neural network. Before applying any model, we have to represent the
sparse system in a form appropriate to input it into a neural network model.

Rose (1972) shows that a sparsity profile of a symmetric matrix can be
described by an undirected graph between $n$ nodes. Starting from this idea,
we define the graph $\mathcal{G}$ associated with a symmetric sparse linear
system.

An undirected graph $\mathcal{G}$ is a pair $(V,E)$ composed of a finite set
of nodes $V$ and a set of undirected edges $E\subseteq\\{\\{u,v\\}\subseteq
V\\}$. A weighted graph additionally associates to every edge $e\in E$ a real
value $w(e)$, called its weight. The adjacency matrix of a weighted undirected
graph is a symmetric matrix $G\in Sym_{n}$ such that
$g_{ij}=g_{ji}=w(\\{i,j\\})$ if $\\{i,j\\}\in E$, otherwise $g_{ij}=0$. Thus,
we can interpret any symmetric sparse matrix as the adjacency matrix of a
weighted undirected graph $\mathcal{G}$.

In our graph representation of a sparse linear system, the coefficients in $A$
become the weights of the undirected edges of the graph $\mathcal{G}$ and
every node of $\mathcal{G}$ represent at the same time a variable of the
system and one of its equations. The order of equations in a linear system
does not influence its solution, so there is no canonical ordering in general.
Nevertheless, when the matrix is symmetric, the order of equations (i.e. the
order of rows in the matrix) and the order of variables (i.e. the order of
columns in the matrix) are linked by the symmetry property of the matrix. To
complete the representation of the system as a graph, we use the elements of b
as node features. In summary, as depicted in Figure 1, the $i^{\text{th}}$
node of the graph $\mathcal{G}$ is described by the feature $b_{i}$ and it is
linked to the $j^{\text{th}}$ node if and only if $a_{ij}\neq 0$. If
$a_{ij}\neq 0$, then the undirected edge $e_{ij}$ is associated to the weight
$a_{ij}$.

Apart from being more suitable to apply deep learning methods to the
resolution of the problem, the graph representation is also inherently
permutation-equivariant. It implies that an identical graph representation
encodes symmetric systems obtained through permutations. This property is
highly desirable because, if $P$ is a permutation matrix, then the system
$A\mathbf{x}=\mathbf{b}$ is equivalent to the system

$(PAP^{\top})\mathbf{y}=P\mathbf{b},\text{ where }\mathbf{y}=P\mathbf{x}\text{
is a permutation of }\textbf{x}.$

On the contrary, the matrix representation of a linear system is not unique,
and the order of rows and columns influences the resolution time of the linear
system. For this reason, permutation heuristics try to find the most efficient
form of a linear system.

Linear systems are also scale-invariant: in Section 4.2, we are going to
explain how to scale edge weights and node features to represent scaled
systems using the same graph.

$\begin{pmatrix}1.0&0.5&0&0&0\\\ 0.5&2.2&4.1&0&1.2\\\ 0&4.1&-1.5&2.0&0\\\
0&0&2.0&3.6&-0.8\\\ 0&1.2&0&-0.8&-0.1\end{pmatrix}\begin{pmatrix}x_{1}\\\
x_{2}\\\ x_{3}\\\ x_{4}\\\
x_{5}\end{pmatrix}=\begin{pmatrix}{\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color<EMAIL_ADDRESS>{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}-1.1}\\\
{\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color<EMAIL_ADDRESS>{\color[rgb]{0.0600000000000001,0.46,1}\definecolor[named]{pgfstrokecolor}{rgb}{0.0600000000000001,0.46,1}\pgfsys@color@cmyk@stroke{0.94}{0.54}{0}{0}\pgfsys@color<EMAIL_ADDRESS>{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}4.8}\end{pmatrix}$

$\big{\downarrow}$

$2.7$$-1.1$$-2.6$$5.4$$4.8$$0.5$$1.0$$4.1$$1.2$$2.2$$2.0$$-1.5$$3.6$$-0.8$$-0.1$

Figure 1: Conversion from the matrix representation to the graph
representation of a sparse symmetric linear system.

##  4 Neural Sparse Solver

Equipped with the graph representation of a sparse linear system, we can now
define our Neural Sparse Linear Solver (NSLS) model. An NSLS is a regression
model that takes as input the graph associated with a sparse linear system
according to the representation introduced in Section 3, and it outputs the
approximate solution of the linear system. Three modules compose an NSLS: a
feature augmentation module, a scaling module, and a graph neural network
module. Such modules interact according to the workflow displayed in Figure 2.

###  4.1 Permutation-Equivariant Feature Augmentation

The naive application of a graph neural network on
$\mathcal{G}_{A,\mathbf{b}}$ is prone to fail because every node stores a
single scalar feature $b_{i}$. This description is insufficient to describe a
node, so we need to enrich the node features and replace the scalar $b_{i}$
with the feature vector $\mathbf{f}_{i}$.

Every permutation-equivariant operation involving $A$ and b that returns a
vector $\mathbf{u}\in\mathbb{R}^{n}$, can be used to augment the features of
the $i^{\text{th}}$ node appending the value $u_{i}$ to it. In practice, all
features are normalized component-wise with the $L_{\infty}$ norm to control
the variance of the input passed to the graph neural network. The $L_{\infty}$
norm is specifically chosen because it is independent of the size $n$.

The most simple augmentation is the diagonal augmentation: the value $a_{ii}$
on the diagonal of the matrix $A$ is added to the feature vector
$\mathbf{f}_{i}$. To design other augmentation operations, we draw inspiration
from existing iterative methods.

Figure 2: Structure of a Neural Sparse Linear Solver. The system is converted
to a graph: the coefficient matrix defines the graph topology, and the
constant term $\mathbf{b}$ gives the node features. Node features are
augmented with permutation-equivariant operations and scaled. A graph neural
network takes this graph as input and predicts the direction of the
approximate solution. Finally, the output of the graph neural network is
scaled to minimize the residual vector.

Jacobi Augmentation. The Jacobi method is an iterative algorithm to solve the
linear system $A\mathbf{x}=\mathbf{b}$ decomposing the matrix as $A=D+M$,
where $D$ is the diagonal matrix, and $M$ is the sum of the lower and upper
triangular parts. Jacobi iteration can be written as

$\begin{cases}\mathbf{x}^{(0)}_{\text{Jacobi}}=\mathbf{0}\\\
\mathbf{x}^{(k+1)}_{\text{Jacobi}}=D^{-1}(\mathbf{b}-M\mathbf{x}^{(k)}_{\text{Jacobi}}).\end{cases}$

We append to the feature matrix $F$ the vectors
$\mathbf{x}^{(1)}_{\text{Jacobi}},\ldots,\mathbf{x}^{(m)}_{\text{Jacobi}}$
after normalization. The normalization here is necessary because the Jacobi
iteration can diverge rapidly if the matrix $A$ is not positive definite or
strictly diagonally-dominant .

We could apply the same idea to similar iterative methods based on other
decompositions of the matrix $A$, like Gauss-Seidel or SOR methods.
Unfortunately, the iterations prescribed by these methods are not permutation-
equivariant because the lower and upper parts of $A$ are different from the
corresponding parts of the permuted matrix $PAP^{\top}$.

Conjugate Gradient Augmentation. A single iteration of the conjugate gradient
method is a permutation-equivariant operation. Therefore, every conjugate
gradient step provides values that can be appended element-wise to the
corresponding node features. We expect this augmentation to be particularly
effective when $A$ is positive-definite, but it is still possible to apply it
to general matrices.

Arnoldi Augmentation. Arnoldi iteration is described by

$\mathbf{x}^{(k)}_{\text{Arnoldi}}=A^{k}\mathbf{b}.$

The GMRES method (Saad & Schultz, 1986) is an iterative method for the
numerical solution of a generic system of linear equations that relies on the
Arnoldi iteration to approximate the solution. Thus, taking inspiration from
GMRES, we enrich the node features appending the normalized version of the
vectors
$\mathbf{x}^{(1)}_{\text{Arnoldi}},\ldots,\mathbf{x}^{(m)}_{\text{Arnoldi}}$.
Here, the permutation equivariance property of the augmentation comes from the
orthogonality of permutation matrices. The normalization is crucial for
Arnoldi augmentation because the repeated multiplication by $A$ can easily
cause overflows or underflows in the values of
$\mathbf{x}^{(k)}_{\text{Arnoldi}}$.

The proposed augmentations can be combined to maximize the accuracy of the
model. The number of elements added to the node features should balance
accuracy gains with the increase in inference time. We compare the end-to-end
performances obtained using every augmentation strategy and their combinations
in Section 5.2.

###  4.2 Scale-Invariant Processing

The graph representation makes the model intrinsically permutation-
equivariant: we now explain how to make the model scale-invariant. The scale
invariance property allows us to describe entire classes of linear systems
with a single representation, and it is even more critical for boosting model
accuracy. Indeed, inputs of uncontrolled scale can hinder the training or make
the model unstable and unreliable. Finally, scale invariance naturally helps
generalization making it possible to manage solution vectors $\mathbf{x}$ of
different scales.

Let us notice that the system $A\mathbf{x}=\mathbf{b}$ is equivalent to

$\mu\frac{A}{||A||}\frac{\mathbf{x}}{||\mathbf{x}||}=\frac{\mathbf{b}}{||\mathbf{b}||},\text{
where }\mu=\frac{||\mathbf{x}||||A||}{||\mathbf{b}||},$

for any norm $||\cdot||$ (or in general for every scaling factor). Let us
denote the scaled quantities above as
$\overline{A},\overline{\mathbf{x}},\overline{\mathbf{b}}$. If the normalized
solution $\overline{\mathbf{x}}$ is known, we can recover

$\mu=\frac{||\overline{\mathbf{b}}||_{2}}{||\overline{A}\overline{\mathbf{x}}||_{2}}.$

If we only know an approximate scaled solution $\widetilde{\mathbf{x}}$, then
we can approximate $\mu$ as

$\widetilde{\mu}=\operatorname*{arg\,min}_{\mu}||\mu\overline{A}\widetilde{\mathbf{x}}-\overline{\mathbf{b}}||_{2}^{2}.$

The minimization problem has the closed-form solution

$\widetilde{\mu}=\frac{\langle\overline{A}\widetilde{\mathbf{x}},\overline{\mathbf{b}}\rangle}{||\overline{A}\widetilde{\mathbf{x}}||_{2}^{2}}.$
(3)

In conclusion, we can represent all scaled versions of a system canonically
using the normalized coefficients $\overline{A}$ and constant terms
$\overline{\mathbf{b}}$. The model receives the system graph with scaled
weights and features and, it predicts a scale-independent approximation of the
solution $\widetilde{\mathbf{x}}$. Hence, we recover the final approximate
solution $\widehat{\mathbf{x}}$ with the formula

$\widehat{\mathbf{x}}=\frac{\widetilde{\mu}||\mathbf{b}||}{||A||}\widetilde{\mathbf{x}}.$
(4)

###  4.3 Graph Neural Network Backbone

The main component of an NSLS is a graph convolutional neural network (GCN)
that maps the scaled augmented node features to the approximate direction of
the solution. The GCN architecture follows the encode-process-decode model
proposed by Battaglia et al. (2018): it is composed of a node encoder,
followed by $m$ residual blocks, and terminates with a solution decoder. The
encoder is a linear layer, while the decoder is 2-layer MLP with LeakyReLU
activation. The decoder outputs a single value per node, and all these values
compose a permutation-equivariant description of the solution vector. For the
residual block, we have tested many layers and their combinations to finally
come up with the structure depicted in Figure 3.

Figure 3: Structure of a residual block in an NSLS.

Residual Connections. Li et al. (2020) empirically find that the pre-
activation variant of residual connections for GCNs, which follows the
ordering Normalization → Activation → GraphConv → Addition, performs better
than the standard ordering introduced by Li et al. (2019). Here, we use a
similar design with an additional GraphConv layer between the normalization
and the activation layers. The supplementary convolutional layer increases the
receptive field of the neural network earlier.

Normalization Layer. We have compared many normalization layers: instance
normalization (Ulyanov et al., 2016), layer normalization (Ba et al., 2016),
graph normalization (Cai et al., 2021), graph size normalization (Dwivedi et
al., 2020), pair normalization (Zhao & Akoglu, 2019). From our experiments, we
see that graph normalization entails consistently the best performances. It is
described by the formula

$\mathbf{x}^{\prime}_{i}=\frac{\mathbf{x}-\alpha\cdot\mathbb{E}[\mathbf{x}]}{\sqrt{\textrm{Var}[\mathbf{x}-\alpha\cdot\mathbb{E}[\mathbf{x}]]+\epsilon}}\cdot\gamma+\beta,$
(5)

where $\alpha,\beta,\gamma$ are scalar parameter learned during training.

Convolution Layer. The most common graph convolution layers, like the ones in
GraphSAGE (Hamilton et al., 2017), GIN (Xu et al., 2018) or GAT (Veličković et
al., 2018), are not suited for weighted graphs. Moreover, the convolutions in
classic graph neural networks like GCN (Kipf & Welling, 2017) or GCNII (Chen
et al., 2020) can lead to zero division errors with some matrices. Among the
graph convolution operators which leverage edge weights, we adopt the
GraphConv layer proposed by Morris et al. (2019). The update formula for the
GraphConv layer is

$\mathbf{x}^{\prime}_{i}=\Theta_{1}\mathbf{x}_{i}+\Theta_{2}\sum_{j\in\mathcal{N}(i)}a_{ji}\mathbf{x}_{j},$
(6)

where $\mathbf{x}_{i}\in\mathbb{R}^{d}$ is the $i$-th node feature vector,
$\Theta_{1},\Theta_{2}\in\mathbb{R}^{d^{\prime}\times d}$ are the parameters
of the layer, and $a_{ji}$ denotes the edge weight from source node $j$ to
target node $i$.

Equation (6) shows that the message aggregation function is the sum. In the
graph neural network literature many aggregation functions are advocated: for
example, mean (Kipf & Welling, 2017), max (Hamilton et al., 2017), softmax and
power mean (Li et al., 2020). Here we choose the sum aggregation function
because, thanks to this choice, the message aggregation reduces to a matrix
multiplication between the scaled coefficient matrix $\overline{A}$ and the
processed node features $F\in\mathbb{R}^{n\times d}$. This choice is
straightforward since the workhorse of many iterative methods is the matrix
multiplication between the coefficient matrix and vectors depending on
$\mathbf{b}$.

###  4.4 Loss Function

The output of the GCN is a scale-independent vector, so we are just interested
in its direction. For this reason we train the model with the cosine distance
loss function:

$\mathcal{L}_{\text{cos}}(\widetilde{\mathbf{x}},\mathbf{x})=1-\frac{\langle\widetilde{\mathbf{x}},\mathbf{x}\rangle}{||\widetilde{\mathbf{x}}||_{2}||\mathbf{x}||_{2}}.$
(7)

Surprisingly, this loss is not sufficient to recover a good approximation. We
have observed that the vectors $\overline{A}\widetilde{\mathbf{x}}$ and
$\overline{b}$ can be near-orthogonal even if the angle between
$\widetilde{\mathbf{x}}$ and $\mathbf{x}$ is tiny. This phenomenon happens
because the transformation induced by an ill-conditioned matrix skews angles.
To counteract this phenomenon, we use a complementary cosine distance loss on
the vectors transformed by $\overline{A}$:

$\mathcal{L}_{\text{res}}(\widetilde{\mathbf{x}},\mathbf{x})=\mathcal{L}_{\text{cos}}(\overline{A}\widetilde{\mathbf{x}},\overline{\mathbf{b}})=1-\frac{\langle\overline{A}\widetilde{\mathbf{x}},\overline{\mathbf{b}}\rangle}{||\overline{A}\widetilde{\mathbf{x}}||_{2}||\overline{\mathbf{b}}||_{2}}.$
(8)

The minimization of $\mathcal{L}_{\text{res}}$ ensures that the projection of
$\overline{A}\widetilde{\mathbf{x}}$ on $\overline{\mathbf{b}}$ does not
collapse into the null vector, making the estimation of $\widetilde{\mu}$
unreliable.

In conclusion, our full loss function is

$\mathcal{L}(\widetilde{\mathbf{x}},\mathbf{x})=\mathcal{L}_{\text{cos}}(\widetilde{\mathbf{x}},\mathbf{x})+\mathcal{L}_{\text{res}}(\widetilde{\mathbf{x}},\mathbf{x}).$
(9)

##  5 Experiments

###  5.1 Results

Dataset. The most used collections of sparse matrices are the Matrix Market
(Boisvert et al., 1997) and the SuiteSparse matrix collection (Davis & Hu,
2011). They contain sparse matrices from many engineering problems, but these
datasets are too small to train a deep learning algorithm on them.
Furthermore, no constant term is reported in conjunction with these matrices.
While simulating the solution is a viable option, it entails unrealistic
linear systems.

The only public dataset of symmetric sparse linear systems we have found is
the StAnD (Grementieri & Finelli, 2022). This dataset contains several solved
linear systems originating in the stability analysis of a structure subjected
to realistic loads. StAnD problems are divided by size: we perform our
experiments on StAnD Small, the subset of linear systems with small
coefficient matrices. Besides, we also conduct our experiments on a synthetic
dataset generated using the same matrices from StAnD Small but with solution
elements uniformly sampled from $[-1,1]$. Our tests reveal very different
behaviors induced by the two datasets, even if they share the same coefficient
matrices. The elements in StAnD Small solution vectors range many orders of
magnitudes (from $10^{-16}$ to $10^{-1}$), and we acknowledge this as the
cause of the discrepancy.

The matrices in StAnD Small have an average dimension of 2115 rows and
columns; this dimension corresponds to the number of nodes in the input system
graph. Whereas the size of the system affects the execution time, the graph
diameter conditions the precision of an NSLS: indeed, every node should
collect information from every other node in its connected component. An NSLS
can gather all values in $\mathbf{b}$, necessary to produce the solution when
the inverse of $A$ is dense, only if this condition holds. Otherwise, if the
diameter of the system graph is larger than the receptive field of the
message-passing neural network at the core of an NSLS, we know that the
network has no theoretical possibility to produce a perfect solution. The
receptive field of graph neural network depends on the connectivity of the
input graph, but in the worst-case scenario of a linear graph it is equivalent
to the number of convolution layers in the GCN. The matrices in StAnD Small
have an average diameter of about 20, so our models should have at least 20
convolution layers to be expressive enough.

Implementation. All models are implemented using PyTorch Geometric (Fey &
Lenssen, 2019). Our models have 10 residual blocks to guarantee a sufficiently
wide receptive field since every residual block contains two GraphConv layers.
The inner feature dimension $d$ is 32 (small) or 128 (medium), and the
training lasts for 50 epochs, or about 315.000 steps using a batch size of 16.
We adopt the Adam optimizer (Kingma & Ba, 2014) with no additional
regularization. Indeed, we find that regularization is detrimental for such a
regression task. We utilize a learning rate $\lambda=10^{-3}$ at the beginning
of the training, dividing it by 10 after 40 and 45 epochs.

Evaluation. We evaluate our models computing on the test set the average of
the mean absolute element-wise error $\epsilon$ and of the relative
approximation error $\delta$ with respect to the ground-truth solution
$\mathbf{x}\in\mathbb{R}^{n}$ defined as

$\epsilon(\widehat{\mathbf{x}},\mathbf{x})=\frac{||\mathbf{x}-\widehat{\mathbf{x}}||_{1}}{n},\qquad\delta(\widehat{\mathbf{x}},\mathbf{x})=\frac{||\mathbf{x}-\widehat{\mathbf{x}}||_{2}}{||\mathbf{x}||_{2}}.$
(10)

The relative error $\delta$ is scale-invariant, allowing us to compare models
trained on different datasets. On the other hand, the mean absolute element-
wise error $\epsilon$ tells us how far the approximation actually is from the
solution. In practical applications, like structural engineering (the domain
of StAnD), we are more interested in the absolute error because the values are
associated with a unit of measure, and tiny errors are negligible.

Table 1 compares our leading models trained on synthetic and StAnD Small data.
The feature augmentations applied to train such models are diagonal and
conjugate gradient augmentations, with 14 steps of the conjugate gradient
method to obtain feature vectors of size 16. This combination of features
attains the best compromise between accuracy and speed. The results on
synthetic data are very satisfactory, exhibiting compelling stability even for
the ill-conditioned matrices typical of structural analysis. The relative
error on actual StAnD data is much higher because the very low norm of the
ground-truth solution distorts it. The absolute error $\epsilon$ tells us
another story, showing that the approximation is adequately good and the NSLS
is sufficiently accurate to replace a classic solver.

Table 1: Evaluation on synthetic and StAnD Small test sets.
Dataset | Inner Size | $\epsilon$ | $\delta$
---|---|---|---
Synthetic | $d=32$ | $9.77\cdot 10^{-3}$ | 3.40%
$d=128$ | $8.03\cdot 10^{-3}$ | 2.80%
StAnD Small | $d=32$ | $8.27\cdot 10^{-6}$ | 18.28%
$d=128$ | $5.60\cdot 10^{-6}$ | 10.55%

###  5.2 Ablation Study

In Table 2, we present the results of the ablation study performed to verify
the effectiveness of each module of an NSLS. To reduce the running time of
experiments, we compare only models with inner feature dimension $d=32$ on
StAnD Small. Even if not specified, we always include the diagonal
augmentation.

Table 2: Ablation study on StAnD Small test set.
Loss | Augmentation | $\epsilon$ | $\delta$
---|---|---|---
$\mathcal{L}_{\textup{cos}}$ | - | $18.3\cdot 10^{-6}$ | 33.93%
$\mathcal{L}_{\textup{cos}}+\mathcal{L}_{\textup{res}}$ | - | $9.01\cdot 10^{-6}$ | 21.66%
$\mathcal{L}_{\textup{cos}}+\mathcal{L}_{\textup{res}}$ | Arnoldi | $8.44\cdot 10^{-6}$ | 19.86 %
$\mathcal{L}_{\textup{cos}}+\mathcal{L}_{\textup{res}}$ | Jacobi | $8.89\cdot 10^{-6}$ | 19.22%
$\mathcal{L}_{\textup{cos}}+\mathcal{L}_{\textup{res}}$ | CG | $\mathbf{8.27\cdot 10^{-6}}$ | 18.28%
$\mathcal{L}_{\textup{cos}}+\mathcal{L}_{\textup{res}}$ | CG + Arnoldi | $8.53\cdot 10^{-6}$ | 19.04%
$\mathcal{L}_{\textup{cos}}+\mathcal{L}_{\textup{res}}$ | CG + Jacobi | $8.30\cdot 10^{-6}$ | 18.12%

The metrics show that the $\mathcal{L}_{\text{res}}$ loss is a fundamental
component of our architecture. Without this term, the model tends to overfit
training data, and the gap between training loss and validation loss increases
during the training.

The comparison between feature augmentation strategies shows that all proposed
methods improve the final result. As expected, the conjugate gradient
augmentation is the most beneficial on positive-definite matrices. On the
other hand, Arnoldi iteration, which does not converge to the solution, is the
least effective strategy, and it is also detrimental when combined with
conjugate gradient.

###  5.3 Comparison with Existing Methods

We compare the running time of NSLS with the fastest open-source C++ libraries
with GPU support for the resolution of linear systems: cuSOLVER (NVIDIA
Corporation, 2022) for direct methods and ViennaCL (Rupp et al., 2016) for
iterative algorithms. Since the matrices in StAnD have the additional property
to be positive-definite, we compare NSLS with algorithms devised especially
for positive-definite systems: Cholesky decomposition (cuSOLVER) and
(Preconditioned) Conjugate Gradient (ViennaCL). We remark that our method does
not need the additional hypothesis of positive-definiteness. On the contrary,
we believe that neural solvers would be much more competitive on general
symmetric systems that are harder to solve with classic algorithms. For
reference, we also measure the execution time of a pure PyTorch implementation
of the conjugate gradient to show that using the same framework and the same
amount of user optimizations, NSLSs are faster.

We combine Cholesky decomposition with the fill-in reducing permutation
strategies SymAMD (Amestoy et al., 1996), SymRCM (George, 1971) or METIS
(Karypis & Kumar, 1998). Similarly, we consider Incomplete Cholesky (Lin &
Moré, 1999) and Incomplete LU (Meijerink & Van Der Vorst, 1977)
preconditioning to foster the convergence of the conjugate gradient method.

Direct methods cannot exchange accuracy for execution time, so they always
reach a very low error. For iterative methods, we record the time needed to
reach a relative error equivalent to the one incurred by the NSLS models
listed in Table 1. The execution time is measured on a GPU NVidia GeForce GTX
1080 Ti.

Table 3: Execution time needed by solvers to reach the average relative error
$\delta$ of the best NSLS with $d=32$ on synthetic and StAnD Small test data.

Algorithm | Variant |  Synthetic Time (ms) |  StAnD Small Time (ms)
---|---|---|---
Cholesky (cuSPARSE) | - | 41.9 | 42.1
SymAMD | 35.3 | 35.3
SymRCM | 36.1 | 36.7
METIS | 31.2 | 31.1
Conjugate Gradient (ViennaCL) | - | 2.16 | 0.84
ICC | 1.97 | 1.20
ILU | 2.49 | 1.77
Conjugate Gradient (PyTorch) | - | 38.1 | 14.4
NSLS | $d=32$ | 13.2 | 13.4
$d=128$ | 18.7 | 18.5

NSLSs are still slower than iterative methods, but they enjoy a better
utilization of GPU resources and a higher degree of parallelization, as can be
seen from the comparison with the pure PyTorch implementation of conjugate
gradient. Moreover, the implementation of NSLS in deep learning frameworks
like PyTorch (Paszke et al., 2019) or TensorFlow (Abadi et al., 2016) makes
these models hardware-agnostic and highly optimized without much effort. We
are sure that the inference time of neural solvers will be further reduced by
the development of technology. Many research papers already propose more
optimized implementations of the sparse matrix multiplication on GPU (Huang et
al., 2020), and deep learning frameworks are actively adding support for
advanced sparse formats and operations. Anyway, many burdens persist: for
example, many neural compilers, like Apache TVM (Chen et al., 2018), do not
support variable-size inputs, and thus the optimization of graph neural
networks is challenging. Furthermore, the ONNX (Bai et al., 2022) support for
sparse operation is very recent (see ONNX opset 16), and accelerators like
ONNX Runtime (ONNX Runtime developers, 2022) still have no support for such
operations.

###  5.4 Integration with Iterative Methods

NSLSs provide a coarse approximation of the solution, so they are suited only
for non-critical problems where the time constraints are more important than
accuracy constraints, for example, in real-time simulations. If more precision
is needed, the output of an NSLS can furnish a high-quality initial guess to a
classic iterative method. Table 4 shows that the NLSL initialization obtained
with the medium model on StAnD Small speeds up the convergence of conjugate
gradient, with a substantial reduction in the number of iterations.

Table 4: The average number of iterations needed by the conjugate gradient
algorithm to reach a certain relative error on StAnd Small using different
initialization vectors.
Relative Error $\delta$ |  Zero Init (iterations) |  NSLS Init (iterations) |  Percentage Reduction
---|---|---|---
1% | 240 | 137 | -43%
0.1% | 362 | 270 | -25%
0.01% | 481 | 395 | -18%
0.001% | 583 | 512 | -12%

##  6 Related Works

Graph Representation. Parter (1961) first introduced explicitly the graph
representation of a sparse matrix. Rose (1972) pursued a detailed graph-
theoretic study of the graph representation for sparse symmetric positive-
definite matrices. The graph representation of the matrix is still exploited
to design permutation algorithms to reduce the number of fill-ins in Gaussian
elimination. This line of work uses an unweighted graph, and nodes have no
features since the graph represents just the matrix of coefficient, not the
entire linear system.

Linear System Solvers. Most of the research on linear system solvers focuses
on reduced classes of linear systems, like symmetric diagonally dominant (SDD)
(Koutis et al., 2012; Spielman & Teng, 2014) or graph-structured (Kyng &
Sachdeva, 2016) linear systems. Recently Peng & Vempala (2021) proposed a
method asymptotically faster than matrix multiplication applicable to the
general case. This last algorithm is comprehensibly very complex, involving
many steps and tricks. Unfortunately, no public implementation of the method
is available, so it is impossible to test it and measure its actual running
time.

Preconditioners. Recent preconditioning methods are specifically designed to
be highly parallel and use efficiently GPUs (Anzt et al., 2018; Dziekonski &
Mrozowski, 2018; Bernaschi et al., 2019; He et al., 2020; Göbel et al., 2021).
Since no preconditioning algorithm works better than others on every matrix,
some authors propose to use deep learning to generate the preconditioner
matrix (Götz & Anzt, 2018; Sappl et al., 2019; Luna et al., 2021).
Preconditioning techniques complement our approach because they improve
existing iterative methods, and we have seen how the integration of NSLSs into
existing algorithms can reduce the number of iterations.

##  7 Conclusion

In this work, we present a deep learning model able to solve approximately
symmetric sparse linear systems. We test the method on sparse positive-
definite linear systems emerging in the field of structural engineering. On
this task, the method is not very accurate and slower than classic methods,
but the results obtained on a synthetic dataset show the potential
capabilities of the model. Moreover, our method suffers from the immaturity of
the support for sparse operations in deep learning frameworks. We are
confident that the execution time of neural solvers will diminish with further
developments of the deep learning libraries.

Overall, our research can be considered a relevant application of a novel
general approach to address numerical analysis problems on sparse matrices.
For example, we could adapt it to predict an optimal preconditioner matrix for
a linear system or an improved initialization for a specific iterative method.
Such examples of further developments show the generality of our approach and
its flexibility. Since the adaptation of algorithms to the sparse case is
challenging and even more troublesome is to parallelize them to use GPUs, we
believe that graph neural networks can be an effective alternative to current
methods.

## References

  * Abadi et al. (2016) Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al.  Tensorflow: a system for large-scale machine learning.  In _Proceedings of the 12th USENIX conference on Operating Systems Design and Implementation_ , pp. 265–283, 2016.
  * Amestoy et al. (1996) Amestoy, P. R., Davis, T. A., and Duff, I. S.  An approximate minimum degree ordering algorithm.  _SIAM Journal on Matrix Analysis and Applications_ , 17(4):886–905, 1996.
  * Anzt et al. (2018) Anzt, H., Huckle, T. K., Bräckle, J., and Dongarra, J.  Incomplete sparse approximate inverses for parallel preconditioning.  _Parallel Computing_ , 71:1–22, 2018.
  * Ba et al. (2016) Ba, J. L., Kiros, J. R., and Hinton, G. E.  Layer normalization.  _arXiv preprint arXiv:1607.06450_ , 2016.
  * Bai et al. (2022) Bai, J., Lu, F., Zhang, K., et al.  Onnx: Open neural network exchange.  https://github.com/onnx/onnx, 2022.  Version: 1.10.2.
  * Battaglia et al. (2018) Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., et al.  Relational inductive biases, deep learning, and graph networks.  _arXiv preprint arXiv:1806.01261_ , 2018.
  * Bernaschi et al. (2019) Bernaschi, M., Carrozzo, M., Franceschini, A., and Janna, C.  A dynamic pattern factored sparse approximate inverse preconditioner on graphics processing units.  _SIAM Journal on Scientific Computing_ , 41(3):C139–C160, 2019.
  * Boisvert et al. (1997) Boisvert, R. F., Pozo, R., Remington, K., Barrett, R. F., and Dongarra, J. J.  Matrix market: a web resource for test matrix collections.  In _Quality of Numerical Software_ , pp. 125–137. Springer, 1997\.
  * Bronstein et al. (2017) Bronstein, M. M., Bruna, J., LeCun, Y., Szlam, A., and Vandergheynst, P.  Geometric deep learning: going beyond euclidean data.  _IEEE Signal Processing Magazine_ , 34(4):18–42, 2017.
  * Cai et al. (2021) Cai, T., Luo, S., Xu, K., He, D., Liu, T.-y., and Wang, L.  Graphnorm: A principled approach to accelerating graph neural network training.  In _International Conference on Machine Learning_ , pp. 1204–1215. PMLR, 2021.
  * Chen et al. (2020) Chen, M., Wei, Z., Huang, Z., Ding, B., and Li, Y.  Simple and deep graph convolutional networks.  In _International Conference on Machine Learning_ , pp. 1725–1735. PMLR, 2020.
  * Chen et al. (2018) Chen, T., Moreau, T., Jiang, Z., Shen, H., Yan, E. Q., Wang, L., Hu, Y., Ceze, L., Guestrin, C., and Krishnamurthy, A.  Tvm: end-to-end optimization stack for deep learning.  _arXiv preprint arXiv:1802.04799_ , 2018.
  * Davies et al. (2021) Davies, A., Veličković, P., Buesing, L., Blackwell, S., Zheng, D., Tomašev, N., Tanburn, R., Battaglia, P., Blundell, C., Juhász, A., et al.  Advancing mathematics by guiding human intuition with ai.  _Nature_ , 600(7887):70–74, 2021.
  * Davis & Hu (2011) Davis, T. A. and Hu, Y.  The university of florida sparse matrix collection.  _ACM Trans. Math. Softw._ , 38(1), 2011.  ISSN 0098-3500.  URL https://doi.org/10.1145/2049662.2049663.
  * Dwivedi et al. (2020) Dwivedi, V. P., Joshi, C. K., Laurent, T., Bengio, Y., and Bresson, X.  Benchmarking graph neural networks.  _arXiv preprint arXiv:2003.00982_ , 2020.
  * Dziekonski & Mrozowski (2018) Dziekonski, A. and Mrozowski, M.  Block conjugate-gradient method with multilevel preconditioning and gpu acceleration for fem problems in electromagnetics.  _IEEE Antennas and Wireless Propagation Letters_ , 17(6):1039–1042, 2018.
  * Fey & Lenssen (2019) Fey, M. and Lenssen, J. E.  Fast graph representation learning with pytorch geometric.  _arXiv preprint arXiv:1903.02428_ , 2019.
  * George (1971) George, J. A.  _Computer implementation of the finite element method_.  Stanford University, 1971.
  * Göbel et al. (2021) Göbel, F., Grützmacher, T., Ribizel, T., and Anzt, H.  Mixed precision incomplete and factorized sparse approximate inverse preconditioning on gpus.  In _European Conference on Parallel Processing_ , pp. 550–564. Springer, 2021.
  * Götz & Anzt (2018) Götz, M. and Anzt, H.  Machine learning-aided numerical linear algebra: Convolutional neural networks for the efficient preconditioner generation.  In _2018 IEEE/ACM 9th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Systems (scalA)_ , pp. 49–56. IEEE, 2018.
  * Grementieri & Finelli (2022) Grementieri, L. and Finelli, F.  Stand: A dataset of linear static analysis problems.  _arXiv preprint arXiv:2201.05356_ , 2022.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, R., and Leskovec, J.  Inductive representation learning on large graphs.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , pp. 1025–1035, 2017.
  * He et al. (2020) He, G., Yin, R., and Gao, J.  An efficient sparse approximate inverse preconditioning algorithm on gpu.  _Concurrency and Computation: Practice and Experience_ , 32(7):e5598, 2020.
  * Huang et al. (2020) Huang, G., Dai, G., Wang, Y., and Yang, H.  Ge-spmm: General-purpose sparse matrix-matrix multiplication on gpus for graph neural networks.  In _SC20: International Conference for High Performance Computing, Networking, Storage and Analysis_ , pp. 1–12. IEEE, 2020.
  * Karypis & Kumar (1998) Karypis, G. and Kumar, V.  A fast and high quality multilevel scheme for partitioning irregular graphs.  _SIAM Journal on scientific Computing_ , 20(1):359–392, 1998.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations (ICLR)_ , 2017.
  * Koutis et al. (2012) Koutis, I., Miller, G. L., and Peng, R.  A fast solver for a class of linear systems.  _Communications of the ACM_ , 55(10):99–107, 2012\.
  * Kyng & Sachdeva (2016) Kyng, R. and Sachdeva, S.  Approximate gaussian elimination for laplacians-fast, sparse, and simple.  In _2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS)_ , pp. 573–582. IEEE, 2016.
  * Lample & Charton (2019) Lample, G. and Charton, F.  Deep learning for symbolic mathematics.  In _International Conference on Learning Representations_ , 2019.
  * Li et al. (2019) Li, G., Muller, M., Thabet, A., and Ghanem, B.  Deepgcns: Can gcns go as deep as cnns?  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 9267–9276, 2019.
  * Li et al. (2020) Li, G., Xiong, C., Thabet, A., and Ghanem, B.  Deepergcn: All you need to train deeper gcns.  _arXiv preprint arXiv:2006.07739_ , 2020.
  * Lin & Moré (1999) Lin, C.-J. and Moré, J. J.  Incomplete cholesky factorizations with limited memory.  _SIAM Journal on Scientific Computing_ , 21(1):24–45, 1999.
  * Logan (2016) Logan, D. L.  _A first course in the finite element method_.  Cengage Learning, 2016.
  * Luna et al. (2021) Luna, K., Klymko, K., and Blaschke, J. P.  Accelerating gmres with deep learning in real-time.  _arXiv preprint arXiv:2103.10975_ , 2021.
  * Meijerink & Van Der Vorst (1977) Meijerink, J. A. and Van Der Vorst, H. A.  An iterative solution method for linear systems of which the coefficient matrix is a symmetric m-matrix.  _Mathematics of computation_ , 31(137):148–162, 1977.
  * Morris et al. (2019) Morris, C., Ritzert, M., Fey, M., Hamilton, W. L., Lenssen, J. E., Rattan, G., and Grohe, M.  Weisfeiler and leman go neural: Higher-order graph neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pp. 4602–4609, 2019.
  * NVIDIA Corporation (2022) NVIDIA Corporation.  cusolver.  https://docs.nvidia.com/cuda/cusolver/, 2022.  Version: 11.6.0.
  * ONNX Runtime developers (2022) ONNX Runtime developers.  Onnx runtime.  https://onnxruntime.ai/, 2022.  Version: 1.10.1.
  * Parter (1961) Parter, S.  The use of linear graphs in gauss elimination.  _SIAM review_ , 3(2):119–130, 1961.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Peng & Vempala (2021) Peng, R. and Vempala, S.  Solving sparse linear systems faster than matrix multiplication.  In _Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA)_ , pp. 504–521. SIAM, 2021.
  * Rose (1972) Rose, D. J.  A graph-theoretic study of the numerical solution of sparse positive definite systems of linear equations.  In _Graph theory and computing_ , pp. 183–217. Elsevier, 1972.
  * Rupp et al. (2016) Rupp, K., Tillet, P., Rudolf, F., Weinbub, J., Morhammer, A., Grasser, T., Jungel, A., and Selberherr, S.  Viennacl—linear algebra library for multi-and many-core architectures.  _SIAM Journal on Scientific Computing_ , 38(5):S412–S439, 2016.
  * Saad & Schultz (1986) Saad, Y. and Schultz, M. H.  Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  _SIAM Journal on scientific and statistical computing_ , 7(3):856–869, 1986.
  * Sappl et al. (2019) Sappl, J., Seiler, L., Harders, M., and Rauch, W.  Deep learning of preconditioners for conjugate gradient solvers in urban water related problems.  _arXiv preprint arXiv:1906.06925_ , 2019.
  * Spielman & Teng (2014) Spielman, D. A. and Teng, S.-H.  Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems.  _SIAM Journal on Matrix Analysis and Applications_ , 35(3):835–885, 2014.
  * Ulyanov et al. (2016) Ulyanov, D., Vedaldi, A., and Lempitsky, V.  Instance normalization: The missing ingredient for fast stylization.  _arXiv preprint arXiv:1607.08022_ , 2016.
  * Van der Vorst (1992) Van der Vorst, H. A.  Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems.  _SIAM Journal on scientific and Statistical Computing_ , 13(2):631–644, 1992.
  * Veličković et al. (2018) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y.  Graph attention networks.  In _International Conference on Learning Representations_ , 2018.
  * Xu et al. (2018) Xu, K., Hu, W., Leskovec, J., and Jegelka, S.  How powerful are graph neural networks?  In _International Conference on Learning Representations_ , 2018.
  * Yannakakis (1981) Yannakakis, M.  Computing the minimum fill-in is np-complete.  _SIAM Journal on Algebraic Discrete Methods_ , 2(1):77–79, 1981.
  * Zhao & Akoglu (2019) Zhao, L. and Akoglu, L.  Pairnorm: Tackling oversmoothing in gnns.  In _International Conference on Learning Representations_ , 2019.