Graph Neural Networks and Applied Linear Algebra
Sparse matrix computations are ubiquitous in scientific computing. With the
recent interest in scientific machine learning, it is natural to ask how sparse
matrix computations can leverage neural networks (NN). Unfortunately,
multi-layer perceptron (MLP) neural networks are typically not natural for
either graph or sparse matrix computations. The issue lies with the fact that
MLPs require fixed-sized inputs while scientific applications generally
generate sparse matrices with arbitrary dimensions and a wide range of nonzero
patterns (or matrix graph vertex interconnections). While convolutional NNs
could possibly address matrix graphs where all vertices have the same number of
nearest neighbors, a more general approach is needed for arbitrary sparse
matrices, e.g. arising from discretized partial differential equations on
unstructured meshes. Graph neural networks (GNNs) are one approach suitable to
sparse matrices. GNNs define aggregation functions (e.g., summations) that
operate on variable size input data to produce data of a fixed output size so
that MLPs can be applied. The goal of this paper is to provide an introduction
to GNNs for a numerical linear algebra audience. Concrete examples are provided
to illustrate how many common linear algebra tasks can be accomplished using
GNNs. We focus on iterative methods that employ computational kernels such as
matrix-vector products, interpolation, relaxation methods, and
strength-of-connection measures. Our GNN examples include cases where
parameters are determined a-priori as well as cases where parameters must be
learned. The intent with this article is to help computational scientists
understand how GNNs can be used to adapt machine learning concepts to
computational tasks associated with sparse matrices. It is hoped that this
understanding will stimulate data-driven extensions of classical sparse linear
algebra tasks.
Nicholas S. Moore, Eric C. Cyr, Peter Ohm, Christopher M. Siefert, Raymond S. Tuminaro
# Graph Neural Networks and Applied Linear Algebra+

###### Abstract

Sparse matrix computations are ubiquitous in scientific computing. Given the recent interest in scientific machine learning, it is natural to ask how sparse matrix computations can leverage neural networks (NN). Unfortunately, multi-layer perceptron (MLP) neural networks are typically not natural for either graph or sparse matrix computations. The issue lies with the fact that MLPs require fixed-sized inputs while scientific applications generally generate sparse matrices with arbitrary dimensions and a wide range of different nonzero patterns (or matrix graph vertex interconnections). While convolutional NNs could possibly address matrix graphs where all vertices have the same number of nearest neighbors, a more general approach is needed for arbitrary sparse matrices, e.g. arising from discretized partial differential equations on unstructured meshes. Graph neural networks (GNNs) are one such approach suitable to sparse matrices. The key idea is to define aggregation functions (e.g., summations) that operate on variable size input data to produce data of a fixed output size so that MLPs can be applied. The goal of this paper is to provide an introduction to GNNs for a numerical linear algebra audience. Concrete GNN examples are provided to illustrate how many common linear algebra tasks can be accomplished using GNNs. We focus on iterative and multigrid methods that employ computational kernels such as matrix-vector products, interpolation, relaxation methods, and strength-of-connection measures. Our GNN examples include cases where parameters are determined _a-priori_ as well as cases where parameters must be learned. The intent with this article is to help computational scientists understand how GNNs can be used to adapt machine learning concepts to computational tasks associated with sparse matrices. It is hoped that this understanding will further stimulate data-driven extensions of classical sparse linear algebra tasks.

os

## 1 Introduction

Artificial intelligence (AI) and machine learning (ML) have drawn a great deal of media attention -- deep fake audio [19], transformer-based AI chatbots [1], stable diffusion AI image generation and even AI Elvis [11] singing a song have all inserted themselves into popular culture. In scientific realms, medical image identification (e.g. detection of cancer) has also generated media attention exposing the potential of ML technologies. But in the realm of applied mathematics the impact of AI/ML has been far less visible to the general public. The goal of this paper is to provide an introduction to graph neural networks (GNNs) and to show how this specific class of AI/ML algorithms can be used to represent (and enhance) traditional algorithms in numerical linear algebra. Associated code which implements the example GNNs is provided at [https://github.com/sandialabs/gnn-applied-linear-algebra/](https://github.com/sandialabs/gnn-applied-linear-algebra/).

Introductions to neural network models in ML and AI typically focus on deep neural networks (DNNs) or convolutional neural networks (CNNs). However, GNNs are can be notably more appropriate for many computational science tasks [29]. DNNs and CNNs generally assume _structured_ input data -- a vector that is a fixed size, or an image where pixels are aligned along Cartesian directions. While structured applications do exist in computational science (uniformly meshed problems for example), many other applications are _unstructured_, again often associated with meshes. Consider a drawing of an object that an engineer wishes to simulate on a computer. After pre-processing, the geometry of the object will be subdivided a mesh. The object may have holes, protrusions or a disparity of features scales such that a structured, Cartesian mesh cannot be generated. As a result, an unstructured mesh is constructed. Application of traditional DNN and CNN tools would likely not be possible for this mesh due to the the unstructured nature of the data. However, GNN models could still be employed on this geometry for associated AI/ML tasks.

So what is a GNN? As suggested by the name, graph neural networks are built on the concept of associating data with edges and vertices of a _graph_. For the above engineered object, the mesh itself can define the graph. Graphs can also represent structured data -- for an image, pixels can be vertices and edges can be used to represent neighboring pixels. The general nature of unstructured graphs make them highly appropriate for a wide range of science and engineering applications. Beyond the meshing example, another example of unstructured data arises from social network graphs where each edge represents a connection between two people. In these type of graph networks, eigenvector calculations provide useful information for determining the influence that a vertex (e.g., a person) has on the rest of the network.
GNNs were first proposed in [16] and [28] as a means of adapting a convolutional neural network to graph problems.

The issue with basic neural networks (e.g. DNNs and CNNs) that GNNs address is that they require a fixed input/feature size. That is, a basic neural network can be viewed as a type of function approximation where the number of input/feature values to the function is always the same when either training the network or when using the network for inference. If a graph has a simple repeatable interaction pattern, then a basic neural network can be effectively applied to a fixed window (i.e., a subset of the graph), which can be moved to address different portions of the network. However, general graphs have no such repeatable pattern. The key difficulty lies in the fact that the number of edges adjacent to each vertex can vary significantly throughout the network. In the case of a social network graph, for instance, a popular person has many friends while a loner might have few interactions. To address this variability, GNNs include the notion of general aggregation functions that allow for a variable number of inputs. Simple aggregation function examples are summation or maximum which are well-defined regardless of the number of inputs. Aggregation functions are used to combine information from edges adjacent to a vertex. These functions can include a fixed number of learnable parameters. For example, an aggregation function might include both a summation and a maximum function whose results are combined in a weighted fashion using a learnable weight. These are then combined with transformation functions associated with edges and vertices, which can also include learnable parameters.

It is well known that many numerical linear algebra algorithms also possess a graph structure. For instance, in Section 3.1 we detail how sparse matrix-vector multiplication transforms matrix entries as graph edge attributes and vector values as graph vertex attributes, into a new set of graph vertex attributes representing the product. This algorithm can be written as a GNN with precise choices of aggregation and transformation functions. We provide additional descriptions of standard iterative methods and components of algebraic multigrid (AMG) preconditioners that can be re-formulated using GNNs. Some of the considered multigrid components include matrix-vector products, Jacobi relaxation, AMG strength-of-connection measures, and AMG interpolation operators. While the associated GNNs do not contain trainable parameters, they illustrate the flexibility of GNNs in representing traditional numerical linear algebra components. The potential for advances in AMG using GNN architectures is made apparent by including trainable parameters in the transformations. These parameters will enable learning complex nonlinear relationships between the feature spaces, implicitly encoding the topological and numerical properties of a sparse matrix. We explore some prototypical use cases in Section 4. We note that more complicated applications of GNNs for multigrid have been considered in the literature and encourage the interested reader to look more deeply [22, 23, 30, 32]. As this paper is educational in nature, we are not proposing new AI/ML powered methods that significantly outperform existing methods, but rather demonstrating how GNNs can allow us to modify existing numerical methods while simultaneously highlighting some of the challenges that must be addressed.

### Developing Intuition: A Viral Example

To make the discussion more accessible, we present an example based on the spread of disease (a subject that we are all unfortunately familiar with). The intent is to help the reader develop intuition about components of the GNN algorithms and data structure. In the text that follows we call out explicitly where we use the example. If the notation and ideas are clear to the reader, these _Viral Interludes_ can be skipped.

Here, we briefly describe the setup of the viral example. We consider a community where each individual interacts on a weekly basis with a fixed set of community members. The one-week period over which the set of interactions occurs is referred to as a cycle. Initially, within the community each individual has a probability of being infected by the disease. Simultaneously, a cure, which also (conveniently) diffuses through individual interactions, has been distributed randomly through the population1. The specific rate of diffusion for the disease and cure, as well as the effectiveness of the intervention is not known in advance, and must be estimated based on observations. The goal then is to develop a model to assess each individual's risk of carrying the infection or carrying the cure after a particular number of cycles using the information about the interactions between individuals.

Footnote 1: Perhaps additional individuals seek treatment after learning about it through peer interactions.

## 2 Graph Neural Network Background

Graph neural networks (GNN) were introduced in [16, 28] to overcome limitations in convolutional neural networks (CNN) that assume a regular structure on the input data to the network. A CNN takes, as input data, structured data that is topologically like a Cartesian mesh, and produces output data with a similarly regular structure. Bitmap images are the usual exemplar for CNN applications. By contrast, the topology of input/output data elements for a GNN is associated with a graph. A GNN takes attributes for each vertex and for each edge and produces a new set of vertex and edge attributes. Additionally, a set of global attributes associated with the entire graph can also be transformed to a new set of global attributes.
Formalizing this description, define a directed graph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\), where \(\mathcal{V}\) is the set of \(N^{v}\) vertices, and \(\mathcal{E}\) is the set of \(N^{e}\) directed edges. A directed edge is defined as an ordered pair of vertices, \((v_{i},v_{j})=e_{ij}\in\mathcal{E}\) where the edge emanates from \(v_{i}\in\mathcal{V}\) and terminates at \(v_{j}\in\mathcal{V}\). The set of features/attributes (data defined on vertices are denoted by \(V=\{\mathbf{v}_{k}\in\mathbb{R}^{n_{e}}:k=1\ldots N^{v}\}\), on edges by \(E=\{\mathbf{e}_{ij}\in\mathbb{R}^{n_{e}}:e_{ij}\in\mathcal{E}\}\), and for the graph are \(\mathbf{g}\in\mathbb{R}^{n_{g}}\). While the terms "features" and "attributes" are both used to describe data associated with edges or vertices in the literature, following Battaglia et al. [4], we will prefer the term attributes. Table 1 summarizes notation used within this paper (including symbols that are introduced shortly).

_Viral Interlude._ For our viral example each vertex in the graph \(v_{i}\in\mathcal{V}\) represents the \(i^{th}\) individual in the community. The initial vertex attribute, prior to the first cycle, is the vector \(\mathbf{v}_{i}\in[0,1]^{2}\in\mathbb{R}^{n_{e}}\) where \(n_{v}=2\). The first component is the probability the individual contains the cure, while the second component is the probability the individual contains the disease. Similarly, the existence of an edge \(e_{ij}\in\mathcal{E}\) indicates an interaction between individual \(i\) and \(j\). Initially, this edge contains as attributes the length of the interaction, which is important in determining both the spread of the infection and the cure, thus \(\mathbf{e}_{ij}\in\mathbb{R}^{n_{e}}\) where \(n_{e}=1\).

### Graph Neural Network Function

A GNN, denoted by \(\mathcal{N}\), is a parameterized function that acts on the attributes of a graph and produces new attributes while maintaining topology:

\[\mathcal{N}\left(\mathcal{G},(V,E,\mathbf{g});\Theta\right)=\{V^{\prime},E^{ \prime},\mathbf{g}^{\prime}\}\,, \tag{1}\]

where \(\Theta\) are the GNN parameters, and a tick denotes the output attributes. Note that the sizes of the attributes can change as a result of the neural network action; for example \(n_{v},n_{e},n_{g}\mapsto n_{v^{\prime}},n_{e^{\prime}},n_{g^{\prime}}\) where \(n_{v}\) can be different from \(n_{v^{\prime}}\), \(n_{e}\) can be different from \(n_{e^{\prime}}\), and \(n_{g}\) can be different from \(n_{g^{\prime}}\). This feature of GNNs will be used in the Viral Interludes, though not in the examples in Section 4. Even if the number of attributes change after the GNN, the topology of the graph where attributes reside is not perturbed by the graph neural network. However, edges and vertices can be marked as on or off in the attributes, though the fundamental topology doesn't change. Also note that the GNN architectures discussed here apply only to directed graphs. Undirected graphs can be represented as well by treating an undirected edge as a pair of directed edges, one in each direction.

In this paper, we focus on GNN models based on message passing that were originally proposed in [13]. This model decomposes the action of the neural network into a sequence of layers, that themselves are

\begin{table}
\begin{tabular}{c l} Symbol & Meaning \\ \hline \(\mathcal{G}\) & Graph \\ \(\mathcal{V}\) & Set of vertices in a graph \\ \(\mathcal{E}\) & Set of edges in a graph \\ \(\mathcal{N}\) & Graph neural network \\ \(\Theta\) & Graph neural network parameters \\ \(N^{v}\) & Number of vertices in a graph: \(|\mathcal{V}|\) \\ \(N^{e}\) & Number of edges in a graph: \(|\mathcal{E}|\) \\ \(V\) & Set of input attributes defined on vertices \\ \(E\) & Set of input attributes defined on edges \\ \(\mathbf{g}\) & Set of input global attributes on a graph \\ \(n_{v}\) & Number of input vertex attributes per vertex \\ \(n_{e}\) & Number of input edge attributes per edge \\ \(n_{g}\) & Number of input global attributes \\ \(v_{k}\) & \(k^{th}\) vertex in the graph \\ \(e_{ij}\) & Directed edge that emanates from \(v_{i}\) and terminates at \(v_{j}\) \\ \(\mathbf{v_{k}}\) & Attributes associated with vertex \(v_{k}\) \\ \(\mathbf{e}_{ij}\) & Attributes associated with edge \(e_{ij}\) \\ \(\overline{\mathbf{e}_{j}}\) & Aggregated attributes from all edges terminating at vertex \(v_{j}\) \\ \(\mathbb{R}^{\text{var}(v)}\) & Variable input space for variable vertex-based functions \\ \(\rho_{e\to v}()\) & Function whose inputs (outputs) are edge (vertex) attributes of a vertex \\ \(\rho_{e\to g}()\) & Function whose inputs (outputs) are edge (global) attributes of the graph \\ \(\rho_{v\to g}()\) & Function whose inputs (outputs) are vertex (global) attributes of the graph \\ \(n_{e\to v}\) & Number of output vertex attributes produced by \(\rho_{e\to v}()\) \\ \(n_{e\to g}\) & Number of output global attributes produced by \(\rho_{e\to g}()\) \\ \(n_{v\to g}\) & Number of output global attributes produced by \(\rho_{v\to g}()\) \\ \end{tabular}
\end{table}
Table 1: Notation summary
parameterized functions of the graph topology and attributes. For instance, a two layer GNN with layers denoted as \(\mathcal{L}_{1}\) and \(\mathcal{L}_{2}\) is

\[\mathcal{N}(\mathcal{G},(V,E,\mathbf{g});\{\Theta_{1},\Theta_{2}\}):=\mathcal{L} _{2}(\mathcal{G},\mathcal{L}_{1}(\mathcal{G},(V,E,\mathbf{g});\Theta_{1}); \Theta_{2}). \tag{2}\]

The specifics of layer architectures are discussed in Section 2.2.

_Viral Interlude._ For our viral model, a GNN function \(\mathcal{N}\) represents the evolution of the disease within this community from the initial state to a final state. \(\mathcal{N}\) is constructed through composition of the layers, where a layer \(\mathcal{L}_{i}\) describes an evolution over a single cycle. In a cycle new attributes for each individual (the set \(V^{\prime}\)) are computed by considering a combination of the vertex attributes \(V\), and the time of interaction specified by the edge attribute \(E\). The parameters \(\Theta\) are not known _a priori_ and must be calibrated based on empirical data about the evolution of the outbreak. (See Section 2.3). In a GNN, the parameter \(\Theta\) encompasses these model unknowns. The layer may also require edge updates based on changing interaction times. Further, if the number of values representing an edge or vertex state changes as a result of applying a layer then we have \(n^{\prime}_{v}\neq n_{v}\) and/or \(n^{\prime}_{e}\neq n_{e}\). For instance, an additional attribute could be added by the GNN layer for each vertex specifying the survival (or not) of an individual. This boolean attribute could have the effect of virtually removing an individual from the graph in subsequent layers of the GNN function.

### Message Passing Layer

While several different variations of message passing GNNs exist, we will describe and utilize the GNN from [4]. The parameterization of a GNN layer is defined by three aggregation functions and three update functions. Because graphs can have varying topology, a method for combining attributes from multiple entities of the same type is necessary. Formally, these aggregation functions are

\[\rho_{e\to v}(\mathbf{e}^{\prime}_{*k}):\mathbb{R}^{\mathrm{var}(v)}\to \mathbb{R}^{n_{e\to v}},\quad\rho_{e\to g}(E^{\prime}):\mathbb{R}^{N^{v}\, \cdot\,n_{e}}\to\mathbb{R}^{n_{e\to g}},\quad\rho_{e\to g}(V^{\prime}): \mathbb{R}^{N^{v}\,\cdot\,n_{v}}\to\mathbb{R}^{n_{e\to g}} \tag{3}\]

where \(\mathbf{e}^{\prime}_{*k}\) denotes the set of updated attributes associated with all edges that terminate at vertex \(v_{k}\). The function \(\rho_{e\to v}()\) is applied at each vertex (i.e., for each \(k=1,\ldots,N_{v}\)) and the notation \(\mathbb{R}^{var(v)}\) denotes that the \(\rho_{e\to v}\) function is variable: the input space is a finite (but variable) number of edge attribute vectors. This function takes the set of attributes from all edges which terminate at a single vertex and aggregates them together into \(n_{e\to v}\) attributes that are collected at this vertex. The other two aggregation functions, \(\rho_{e\to g}\) and \(\rho_{v\to g}\), are similar but combine updated attributes \(E^{\prime}\) associated with all graph edges and updated attributes \(V^{\prime}\) associated with all graph vertices respectively. Some simple examples of aggregation functions include summation, minimum, and maximum. The aggregated attributes are then used to update the attributes of the vertices and the graph.

The update functions are

\[\phi_{e}(\mathbf{e}_{ij},\mathbf{v}_{i},\mathbf{v}_{j},\mathbf{g }):\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{v}}\times\mathbb{R}^{n_{e}}\times \mathbb{R}^{n_{g}}\to\mathbb{R}^{n^{\prime}_{e}}, \tag{4}\] \[\phi_{v}(\mathbf{v}_{i},\overline{\mathbf{e}}_{i},\mathbf{g}): \mathbb{R}^{n_{v}}\times\mathbb{R}^{n_{e\to v}}\times\mathbb{R}^{n_{g}} \to\mathbb{R}^{n^{\prime}_{e}},\] (5) \[\phi_{g}(\mathbf{g},\rho_{e\to g}(E^{\prime}),\rho_{v\to g}(V^{ \prime})):\mathbb{R}^{n_{g}}\times\mathbb{R}^{n_{e\to g}}\times \mathbb{R}^{n_{v\to g}}\to\mathbb{R}^{n^{\prime}_{g}} \tag{6}\]

where \(\overline{\mathbf{e}}_{i}\) are edge attributes that have been aggregated to vertex \(i\) (i.e., \(\overline{\mathbf{e}}_{i}=\rho_{e\to v}(\mathbf{e}^{\prime}_{i})\)). Each of these functions (with any implied parameters), takes, as input, the attributes of a graph entity and the attributes of entities in its neighborhood and transforms them into an updated attribute associated with the original graph entity. For instance \(\phi_{e}\) is associated with updating an edge's attributes in the graph. The edge's neighborhood includes two vertices, and the graph itself. Therefore, it takes as input the attributes of the edge, its two neighboring vertices and the graph's attributes. The vertex update, \(\phi_{v}\), is associated with a vertex in the graph. The neighborhood of the vertex includes the graph itself and all the edges which terminate at the vertex. The number of terminating edges can vary from vertex to vertex, which necessitates use of the variadic aggregation function \(\rho_{e\to v}\) that combines all the edge attributes together. The original attributes of the vertex, the aggregated edge attributes, and the global attributes are then used as arguments to the vertex update function. A similar process occurs in the global update \(\phi_{g}\) as well, where all updated edge attributes are aggregated together and all updated vertex attributes are aggregated together and used to update the global graph attributes.

With the the aggregation and update functions described, we can now explain how they are combined in Algorithm 1 to compute the action of a GNN layer. The first **for** loop in the algorithm transforms input edge attributes using the update function \(\phi_{e}\). This is depicted graphically in Figure 1 with the initial state of the graph shown on the left graph and the final state shown on the right. For a single edge \(e_{12}\), its attributes, along with those of the neighboring vertices, and the graph are the input to \(\phi_{e}\) yielding an updated edge attribute \(\mathbf{e}^{\prime}_{12}\). While only one update is shown, all edges are updated (potentially in parallel) using the same update function.
```
0: Graph \(G\), vertex attributes \(V=\{\mathbf{v}_{j}:j=1,\ldots,N^{v}\}\), edge attributes \(E=\{\mathbf{e}_{ij}:e_{ij}\in\mathcal{E}\}\) and global attributes \(\mathbf{g}\) Output: updated edge attributes \(E^{\prime}\), vertex attributes \(V^{\prime}\) and global attributes \(\mathbf{g}^{\prime}\)
```

**Algorithm 1** Computation of a Graph Network Layer

With the edges transformed, the next **for** loop in Algorithm 1 updates the vertex attributes. Each iteration contains two steps, the first being the aggregation of all edges terminating at the current vertex using \(\rho_{e\to v}\), and the second updating the vertex's attribute. This is depicted in Figure 2 for \(v_{4}\). Here, the edge neighborhood contains edges that terminate at \(v_{4}\). The attributes \(\mathbf{e}_{24}^{\prime},\mathbf{e}_{34}^{\prime}\) and \(\mathbf{e}_{14}^{\prime}\) for these three edges are aggregated using \(\rho_{e\to v}\) yielding \(\bar{\mathbf{e}}_{4}\). The aggregated edge attributes, with the vertex attribute \(\mathbf{v}_{4}\), and the global attribute \(\mathbf{g}\) are input to the update function \(\phi_{v}\). While only one update is shown, this procedure is repeated (potentially in parallel) until all vertex attributes are updated.

The final phase of the GNN layer is to transform the global attributes of the graph. This is done by first aggregating all the transformed edge attributes using \(\rho_{e\to g}\) and all the transformed vertex attributes using \(\rho_{v\to g}\). Then the new global attributes \(\mathbf{g}^{\prime}\) are computed using \(\phi_{g}\). The transformed edge, vertex, and graph attributes are returned in the final line of Algorithm 1.

A notable quality of the message passing layer is that if the graph attributes are ignored (\(\mathbf{g}\) doesn't change), then the action of the layer only considers the one-ring (or distance-one) neighborhood of a vertex for attribute update. So the action of two consecutive GNN layers will update attributes based on vertices two edges away, or the two-ring neighborhood. This breadth-first approach has implications on the use of GNNs for processing subgraph information independently, leaving intriguing possibilities for parallelism. Additionally, each pass of the layer results in an exchange of information over attributes within a single vertex one-ring neighborhood. Consequently, multiple applications of the layer communicate to ever broader vertex neighborhoods.

Figure 1: Edge feature update for \(\mathbf{e}_{12}\). The concatenated attributes for the edge, neighboring vertices, and global graph are input to the vertex update function \(\phi_{e}\). The output \(\mathbf{e}_{12}^{\prime}\) corresponds to updated edge attributes. This process (not shown) is repeated for all edges in parallel, using the same update function.
_Viral Interlude._ The update for the viral outbreak graph described conceptually above can be made (more) formal by defining appropriate aggregation and update functions. The edge update function \(\phi_{e}\) must provide a measurement of an interaction contributing to the cure and disease state of an individual. Further, if a graph function is to be applied multiple times then the time of interaction must also be represented in the results of this function:

\[\mathbf{e}^{\prime}_{ij}\leftarrow\phi_{e}(\mathbf{e}_{ij},\mathbf{v}_{i}, \mathbf{v}_{j},\mathbf{g})=[\mathbf{e}_{ij},\mathbf{i}_{ij}] \tag{7}\]

where \(\mathbf{e}_{ij}\) is the original attribute, the time of interaction, and \(\mathbf{i}_{ij}\) is a quantification of the impact of the interaction (note that in the context of neural networks the quantity \(\mathbf{i}_{ij}\) is not understood in a precise way, rather the GNN will _learn_ to interpret it). For this layer \(n^{\prime}_{e}=n_{e}+1\). While this quantifies the impact of a single interaction, the end state of an individual is based on all their interactions specified by connected graph edges. This is the intent of aggregation function, \(\rho_{e\to v}\). For example, the sum of the interaction risks might be a useful metric, as would the maximum and/or minimum risk, depending on the dynamics of the cure and disease. More formally, a vector-valued function that quantifies the external risk to an individual \(j\) over all interactions is

\[\bar{\mathbf{e}}_{j}=\rho_{e\to v}(\mathbf{e}^{\prime}_{ij})=\left[\sum_{i} \mathbf{e}^{\prime}_{ij}\ \,\max_{i}\mathbf{e}^{\prime}_{ij},\ \ \min_{i}\mathbf{e}^{\prime}_{ij}\right]. \tag{8}\]

With the external risks to an individual quantified, the vertex attributes specifying the cure/disease probability for the individual must be updated using \(\phi_{v}\).

### Training

Typically in GNN applications, MLP neural networks are used as the update functions. Each of these MLPs is parameterized by a set of weights and biases. Therefore, the parameters of a GNN is the union of the weights and biases of all the update functions. We denote the collected set of parameters for the neural network as \(\Theta\). The activity of tuning these parameters to learn to perform a set task is referred to as training. The are different types of tasks that can be learned, but the most common is referred to as supervised learning. In this case, a desired mapping is specified empirically using a dataset composed of input-output pairs. For exposition, in this section we define the data set as \((X_{k},Y_{k})\) for \(k=1\ldots n_{d}\), where \(X_{k}\) is the domain element mapping to \(Y_{k}\), the range element. For GNN's, the domain and range values are defined as vertex, edge and/or global attributes. The parameters of a neural network are said to be trained when a _loss_ function is minimized. An idealized supervised learning problem and loss function is

\[\Theta^{*}=\operatorname*{arg\,min}_{\Theta}\left(\mathcal{L}(\Theta):=\frac{ 1}{n_{d}}\sum_{k=1}^{n_{d}}l_{\Theta}(\mathcal{N}(X_{k};\Theta),Y_{k})\right) \tag{9}\]

where the function \(l_{\Theta}\) measures the difference between the GNN prediction for the attributes \(X_{k}\) and the empirical target \(Y_{k}\). Choices for \(l_{\Theta}\) include norms for regression, or cross entropy for classification.

With the loss function selected, an algorithmic approach to find the optimal parameters \(\Theta^{*}\) is required. Initially the weights and biases are selected at random (see [14, 15, 8, 17] for example approaches and considerations). Then an iterative method is used to incrementally reduce the loss function and improve

Figure 2: Vertex feature update for \(\mathbf{v}_{4}\). The updated attributes for incoming edges are aggregated using \(\rho_{e\to v}\) yielding aggregated edge information \(\bar{\mathbf{e}}_{4}\). The aggregated attributes are concatenated with vertex attributes and global attributes to define the input to the vertex update function \(\phi_{v}\). The output \(\mathbf{v}_{4}^{\prime}\) corresponds to updated vertex attributes. This process is repeated for all vertices in parallel (not shown), using the same aggregation and update functions.
the prediction by adjusting the GNN parameters. The most common iterative optimization methods for training neural networks are gradient descent algorithms. At each iteration, the algorithm calculates the loss for the current \(\Theta\) parameters by forward propagation through a GNN and evaluating a loss on the predicted versus observed output values. In a second step of the iteration, the gradient of the loss with respect to \(\Theta\) (e.g. \(\nabla_{\Theta}\mathcal{L}\)) is computed using the celebrated _backpropagation_ algorithm [15]. Simplistically, the negative of the gradient is used to update the parameters and the iteration is repeated. An important metric to measure the rate of the reduction of the loss is the _epoch_. One epoch corresponds to each entry in the dataset having been used in a gradient descent step once. Values of tens to thousands of epochs are not uncommon. While gradient descent [24] is a common algorithm in training, more sophisticated algorithms can lead to more robust results, for example Stochastic Gradient Descent [26], RMSProp [33], and Adam [20]. For a review of optimization methods for machine learning see [5, 15].

The ability to calculate the gradient for a large number of parameters is automated by software libraries which perform automatic differentiation. Common libraries for machine learning such as PyTorch [25], TensorFlow [2], Jax [6], and others all include this capability, allowing gradient descent algorithms to be applied without requiring user-specified gradient functions.

_Viral Interlude._ Applied to our viral example, the goal of training is to determine parameters that define the transmission rates, and effectiveness of the cure versus the spread of the disease. In this supervised learning example, the dataset input is defined as the initial distribution of the virus and cure over the vertices of the graph, while the output data will be the distribution of the virus and cure after multiple cycles. The dataset would be comprised of graphs defined by multiple communities already inflicted by the virus (thus an observation can be obtained). The model edge update function is parameterized by a neural network that defines the transfer of cure/disease through an individual interaction. The vertex update is also a neural network model that defines the uptake of the cure/disease by a single individual through multiple interactions. Prior to training, the neural network contains only the context of the interactions, not the specifics of the cure/viral diffusion. These specifics would be learned by defining a regression-based loss on the vertex attributes associated with the probability of an individual having the virus/cure. Once trained, the GNN model can be applied to communities (graphs) that have an initial viral load to predict a final state distribution.

## 3 Linear Operations as Graph Network Layers

This section demonstrates how linear algebra computations can be carried out utilizing the language and structure of graph neural networks. These familiar algorithms are selected to make the graph neural network framework described above concrete. We note that the representation of these operations as graph neural networks are not unique, and multiple representations may exist for a given calculation. As such, this is not meant to be an exhaustive list, but instead a general guide to setting up and representing different linear algebra operations as graph neural networks.

We begin by briefly illustrating the natural relationship between sparse matrices and graphs. A weighted graph with \(n\) vertices can represent a square \(n\times n\) matrix \(A\). Each nonzero \(A_{ij}\) defines a weight for a directed edge that emanates from vertex \(j\) and terminates at vertex \(i\). Notice that this definition uses self edges to represent entries on the matrix diagonal. Alternatively, one can omit self edges and instead store \(A_{ii}\) at the \(i^{th}\) vertex. The inclusion or omission of self edges affects how linear algebra operations are represented and how data propagates through the network as different information is available during specific update/aggregation phases. The algorithms in this section assume that the graph neural network includes self edges, except for the matrix-vector example in Section 3.1. Additionally, our edge orientation choice is natural for matrix-vector products \(y=Ax\) where the \(i^{th}\) entry of \(y\) is defined by

\[y_{i}=\sum_{j\in A_{i*}}A_{ij}x_{j}\]

and \(A_{i*}\) denotes the set of all nonzero entries in the \(i^{th}\) matrix row. Here, information from neighboring edges must be gathered at the \(i^{th}\) vertex. This is accomplished by an aggregation function as the edges in the \(i^{th}\) matrix row terminate at the \(i^{th}\) vertex. More generally, our edge orientation choice emphasizes that information that directly influences vertex \(i\) flows into vertex \(i\). Finally, we note that we only consider square matrices. However, non-square matrices can be represented as bipartite graphs where rows and columns each have their own distinct set of vertices. Bipartite graphs extensions of Algorithm 1 are possible, but are beyond the scope of this educational survey.

We start with foundational linear algebra algorithms: sparse matrix-vector product and a matrix-weighted norm. Next, we describe three simple iterative methods: a weighted Jacobi linear solver, a Chebyshev linear solver, and a power method eigensolver. We conclude with some kernels used within an advanced algebraic multigrid linear solver. All of these examples follow Algorithm 1 for the structure of the graph neural network. The details of each operation specific component, such as the update functions (\(\phi_{v}\), \(\phi_{e}\), \(\phi_{g}\)) and aggregation functions (\(\rho_{e\to v}\), \(\rho_{e\to g}\), \(\rho_{v\to g}\)), are summarized in Tables 1-1.
_Code Availability_. Code for each of the GNNs given in this chapter can be found in the repository located at [https://github.com/sandialabs/gmn-applied-linear-algebra/](https://github.com/sandialabs/gmn-applied-linear-algebra/). Code is provided both in MATLAB script and python via PyTorch and the PyTorch Geometric package. Each example provides the code implementation of the layer as well as a small demonstration which shows the output of the GNN matches that of the "traditional" method.

### Sparse Matrix-Vector Product

Sparse matrix-vector products are fundamental building blocks of many numerical linear algebra algorithms, such as Krylov methods. Table 1 illustrates a GNN that computes \(y=Ax\).

The left block of the table describes the input, intermediate and output data for each graph component (vertex, edge, global). The right column is subdivided into the update functions required by Algorithm 1 on the left, and the aggregation functions. Further, if there are distinct layers in the network, the right column lists the sets of update and aggregation functions.

For the sparse matrix vector product the nonzero entries \(A_{ij}\) naturally correspond to edges, so we assign them to edges as fixed (i.e., unmodified by the GNN) data objects. The vector entries \(x_{i}\) and \(y_{i}\) are assigned to vertices, with \(x_{i}\) fixed and \(y_{i}\) mutable (i.e., modified by the GNN) and initialized to zero. The upper left graph in Figure 1 illustrates the starting condition of the GNN with data placed according to the description in Table 1.

We now follow the update and aggregation functions of Algorithm 1 in order. Each step is a transition between graph states in Figure 1. The edge feature update \(\phi_{e}\) takes fixed edge feature \(A_{ij}\) and fixed vertex feature \(x_{j}\) and multiplies them, storing the result on each edge. Summation is employed for edge aggregation, denoted \(\phi_{e}\), with the result stored on the vertex in the vertex update phase, \(\phi_{v}\). The lower left graph in Figure 1 shows the final state of the GNN, with the entries \(y_{i}\) containing the matrix vector product located at graph vertices. In Table 1 and all the tables that follow, it is understood that operations such as \(A_{ij}x_{j}\) and \(\sum_{j}c_{ij}\) only operate on nonzero entries of \(A\) and \(c\) respectively. That is, the edges must be included in the graph. The table can be modified if self-edges are not stored. To that end, Table 2 shows a slightly more complicated GNN. The difference is that \(A_{ii}\) is stored on vertices and so the \(A_{ii}x_{i}\) term is added in the vertex update function \(\phi_{v}\).

### Matrix-Weighted Norm

Similar to the matrix-vector product, a GNN for computing a matrix-weighted vector norm, \(\left\|x\right\|_{W}=\sqrt{x^{T}Wx}\), is shown in Table 3.

The matrix-vector product from Table 1 is modified in the vertex update and adds operations for vertex-to-global aggregation and global

\begin{table}
\begin{tabular}{|l l l||l l|l l|} \hline \multicolumn{3}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\ \multicolumn{3}{|c||}{Fixed} & \multicolumn{1}{l||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{2}{c|}{Aggregation} \\ \hline \hline \multirow{3}{*}{Edge} & \(A_{ij}\) & \(c_{ij}\) & \(\phi_{e}\) & \(c_{ij}=A_{ij}x_{j}\) & \(\rho_{e\to v}\) & \(\overline{c}_{i}=\sum_{j}c_{ij}\) \\  & \(x_{i}\) & \(y_{i}\) [Output] & \(\phi_{v}\) & \(y_{i}=\overline{c}_{i}\) & \(\rho_{e\to g}\) & — \\ \cline{1-1}  & \(\quad\) — & \(\quad\) — & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 1: Sparse Matrix-Vector Product as a Graph Network with Self-Edges

\begin{table}
\begin{tabular}{|l l l||l l|l l|} \hline \multicolumn{3}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\ \multicolumn{3}{|c||}{Fixed} & \multicolumn{1}{l||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{2}{c|}{Aggregation} \\ \hline \hline \multirow{3}{*}{Edge} & \(A_{ij}\) & \(c_{ij}\) & \(\phi_{e}\) & \(c_{ij}=A_{ij}x_{j}\) & \(\rho_{e\to v}\) & \(\overline{c}_{i}=\sum_{j}c_{ij}\) \\  & \(x_{i}\) & \(y_{i}\) [Output] & \(\phi_{v}\) & \(y_{i}=\overline{c}_{i}\) & \(\rho_{e\to g}\) & — \\ \cline{1-1}  & \(\quad\) — & \(\quad\) — & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 2: Sparse Matrix-Vector Product as a Graph Network without Self-Edges
update. After the edge update and the edge-to-vertex aggregation, \(\overline{e}\) contains the vector \(Wx\). Therefore, in the vertex update, we multiply by \(x_{i}\) to obtain the vector \(y_{i}\). In the vertex-to-global aggregation, \(y_{i}\) is summed to obtain \(x^{T}Wx\). Finally, in the global update, the square root is taken to yield the vector norm.

### Weighted Jacobi Iteration

Weighted Jacobi iteration is a simple method to solve a linear system \(Ax=b\) for \(x\). One iteration of weighted Jacobi is written as an update formula

\[x^{k+1}=x^{k}+\omega D^{-1}(b-Ax^{k}), \tag{10}\]

where \(k\) is the iteration index, \(D\) is the matrix diagonal, \(b\) is the right-hand side, \(x^{k}\) is the solution at the \(k^{th}\) Jacobi iteration, and \(\omega\) is the weight parameter. For the GNN representation of the weighted Jacobi method, the edge data includes the matrix entries \(A_{ij}\). The initial vertex data is comprised of the matrix diagonal, \(A_{ii}\), and the right-hand side vector, \(b_{i}\). The weight parameter, \(\omega\) is included as fixed global data. Note that while the matrix diagonal, \(A_{ii}\), is stored as vertex data, this GNN layer uses self-edges as well, so the matrix diagonal is also accessible as regular edge data. This is done to allow easy access to this data for both the matrix-vector product and the Jacobi update. For the functions, all of the components of the sparse matrix-vector product are retained, except for the vector update, \(\phi_{v}\). That is replaced to use the aggregated edge data inside a Jacobi-style update, (10) to yield the GNN layer shown in Table 4.

### Chebyshev Iterative Solver

The Chebyshev method is another simple iterative scheme. It can be viewed as a type of multi-step Jacobi algorithm where different weights are used for each step and these weights are optimal in some sense given knowledge about the largest and smallest eigenvalue (\(\lambda_{max}\) and \(\lambda_{min}\)) of the symmetric positive definite matrix \(A\). An example algorithm is given in Algorithm 2, following the algorithm presentation in [27].

```
1:\(\mathcal{L}=\{x^{k},y^{k}\}\), \(\overline{
For the Chebyshev GNN, most of the computation is nodal and the GNN primarily serves to propagate residual information. When multiple values are updated on a single line of an update function, the values update from left to right, so parameters are updated in the correct order. In each iteration of the Chebyshev algorithm the GNN uses Layer 1 to update the current iterate. Then Layer 2 updates the residual information. Finally, Layer 3 updates the direction vector based on the the updated residual. If \(N\) iterations of the Chebyshev method are desired, run the GNN in Table 3.5\(N\) times consecutively.

### Power Method

Since the power method consists for symmetric matrices consists of matrix-vector multiplications and norms, we can represent the power method as a GNN as well. The details are given in Tables 3.6 and 3.7. We split the algorithm into two networks. The network in Table 3.6 is the iterative portion which updates the \(b\) vector to an improved approximation of the eigenvector associated with the dominant eigenvalue. The second network computes the Rayleigh quotient to give the eigenvalue estimate. Thus, for more than a single iteration of the power method, the network in Table 3.6 is applied successively for the desired number of iterations, and then the Rayleigh quotient is computed from the network in Table 3.7.

\begin{table}
\begin{tabular}{|l l l||l l|l l|} \hline \multicolumn{4}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\ \multicolumn{2}{|c|}{Fixed} & \multicolumn{2}{c||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{2}{c|}{Aggregation} \\ \hline \hline \multicolumn{4}{|c||}{Edge} & \multicolumn{1}{c}{\(A_{ij}\)} & \multicolumn{1}{c}{\(c_{ij}\)} & \multicolumn{1}{c}{\(\phi_{e}\)} & \multicolumn{1}{c|}{\(c_{ij}=A_{ij}x_{j}\)} & \multicolumn{1}{c|}{\(\rho_{e\to v}\)} & \multicolumn{1}{c|}{\(\overline{c}_{i}=\sum_{j}c_{ij}\)} \\ \multicolumn{2}{|c|}{Vertex} & \multicolumn{1}{c}{\(A_{ii},b_{i}\)} & \multicolumn{1}{c}{\(x_{i}\) [Output]} & \(\phi_{v}\) & \(x_{i}=x_{i}+\omega A_{ii}^{-1}(b_{i}-\overline{c}_{i})\) & \(\rho_{e\to g}\) & — \\ \multicolumn{2}{|c|}{Global} & \multicolumn{1}{c}{\(\omega\)} & — & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 3.4: Weighted Jacobi Iteration as a Graph Network

\begin{table}
\begin{tabular}{l l l} \hline \multicolumn{2}{|c||}{**function** (\(x\)) = Cheby(\(A,b,x_{0},\lambda_{max},\lambda_{min}\),\(N\))} \\ \multicolumn{2}{|c|}{\(x\gets x_{0}\) ; } & \multicolumn{1}{c}{\(r\gets b-Ax\)} \\ \multicolumn{2}{|c|}{\(\theta\leftarrow(\lambda_{max}+\lambda_{min})/2\) ; } & \multicolumn{1}{c}{\(\delta\leftarrow(\lambda_{max}-\lambda_{min})/2\) ; } & \multicolumn{1}{c}{\(\sigma\leftarrow\theta/\delta\)} \\ \multicolumn{2}{|c|}{\(\rho\leftarrow 1/\sigma\) ; } & \multicolumn{1}{c}{\(d\leftarrow(1/\theta)\)\(r\)} \\ \multicolumn{2}{|c|}{**for**\(i\in\{1\ldots N\}\)**do**} \\ \multicolumn{2}{|c|}{\(x\gets x+d\)} \\ \multicolumn{2}{|c|}{\(r\gets r-Ad\)} \\ \multicolumn{2}{|c|}{\(\rho_{\text{prior}}\leftarrow\rho\)} \\ \multicolumn{2}{|c|}{\(\rho\gets 1/(2\sigma-\rho)\)} \\ \multicolumn{2}{|c|}{\(d\leftarrow\rho\)\(\rho_{\text{prior}}\)\(d\)\(+\)\((2\rho/\delta)\)\(r\)} \\ \multicolumn{2}{|c|}{**end for**} \\ \multicolumn{2}{|c|}{**return**\(x\)} \\ \multicolumn{2}{|c|}{**end function**} \\ \hline \multicolumn{4}{|c||}{_Chebyshev Solver Iteration as a Graph Network_} \\ \hline \multicolumn{4}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\ \multicolumn{2}{|c|}{Fixed} & \multicolumn{1}{c||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{1}{c|}{Aggregation} \\ \hline \hline \multicolumn{4}{|c||}{Edge} & \multicolumn{1}{c}{\(A_{ij}\)} & \multicolumn{1}{c}{\(c_{ij}\)} & \multicolumn{1}{c}{\(\phi_{e}\)} & — & \(\rho_{e\to v}\) & — \\ \multicolumn{2}{|c|}{Vertex} & \multicolumn{1}{c}{\(r_{i},d_{i},x_{i}\) [Output]} & \(\phi_{v}\) & \(x_{i}=x_{i}+d_{i}\) & \(\rho_{e\to g}\) & — \\ \multicolumn{2}{|c|}{Global} & \multicolumn{1}{c}{\(\delta\), \(\sigma\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} \\ \multicolumn{2}{|c|}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} \\ \multicolumn{2}{|c|}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} \\ \multicolumn{2}{|c|}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} \\ \multicolumn{2}{|c|}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior}}\)} & \multicolumn{1}{c}{\(\rho,\rho_{\text{prior
### Advanced Algebraic Mutlgrid Solver

We conclude this section with two algebraic multigrid kernels: one for interpolation and another associated with coarsening a matrix graph. Algebraic multigrid (AMG) is an advanced linear solver. In the description that follows, we consider solving the matrix system \(A_{1}x_{1}=b_{1}\) for the unknown vector \(x_{1}\). While Jacobi and Chebyshev may require a very large number of iterations to reach an acceptable solution, AMG can converge rapidly on many important matrix systems. For example, AMG is known to converge efficiently on matrices arising from discrete approximations of elliptic partial differential equations (PDEs) such as the heat equation [7]. In particular, the number of required iterations can be independent of the dimension of the matrix \(A_{1}\). Thus, increasingly larger matrix systems do not necessarily require increasingly more iterations. As matrix systems can include over a billion unknowns, this property is extremely advantageous. AMG methods are primarily used to solve discrete versions of PDEs (including many non-elliptic PDEs), though they have been applied to other systems such as those arising from certain network problems. While a full understanding of AMG is beyond our scope (see [7] for a more thorough explanation of the method), we present a few AMG ideas using a simple example.

Consider a metal plate with a hole where the plate is clamped at two locations. One clamp is fixed at \(40\lx@math@degree C\) while the other clamp is fixed at \(10\lx@math@degree C\). Suppose a discrete representation of Poisson's equation is used to model the steady-state heat distribution within the plate where the associated matrix \(A_{1}\) only has nonzeros on the diagonal and for each edge of a mesh. The vector \(b_{1}\) contains only a few nonzeros to represent the fixed temperature of the plate where it is held by the clamps. The other boundary conditions indicate that the heat gradient normal to the plate boundaries is zero. When starting with an initial guess of \(x_{1}^{(0)}=0\), the first Jacobi iteration produces \(x_{1}^{(1)}=\omega D_{1}^{-1}A_{1}b_{1}\). Due to the sparsity pattern of \(A_{1}\), \(x_{1}^{(1)}\) must only have a few nonzeros, all within a graph distance of one from one of the two clamp locations. In general, the \(k^{th}\) Jacobi iteration extends \(x_{1}^{(k-1)}\)'s nonzero regions along edges adjacent to these regions. Thus, one can see that many iterations are required for the influence of the two clamps to propagate throughout the entire plate, and so it should not be surprising that many further iterations are required to reach a converged solution. As the mesh is refined, more iterations are needed to propagate information from the clamps to the rest of the plate. The situation is similar for \(N\)-step Chebyshev solvers, which propagate information \(N\) times faster than Jacobi but are also \(N\) times as

\begin{table}
\begin{tabular}{|l l l||l||l|l|} \hline \multicolumn{3}{|c||}{} & \multicolumn{2}{c||}{Data} & \multicolumn{2}{c|}{Functions} \\ \multicolumn{3}{|c||}{Fixed} & \multicolumn{1}{c||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{2}{c|}{Aggregation} \\ \hline \hline \multirow{3}{*}{Edge} & \(A_{ij}\) & \(c_{ij}\) & \(\phi_{e}\) & \(c_{ij}=A_{ij}b_{j}\) & \(\rho_{e\to v}\) & \(\overline{c}_{i}=\sum_{j}c_{ij}\) \\  & — & \(b_{i},y_{i}\) & \(\phi_{v}\) & \(b_{i}=\overline{c}_{i}\) & \(\rho_{e\to g}\) & — \\  & — & \(h_{1},n_{A},\lambda_{\max}[Output]\) & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \cline{3-6}  & & & & \multicolumn{2}{c|}{Layer 2 (\(h_{1}=\|b\|\))} \\ \cline{3-6}  & & & & \multicolumn{2}{c|}{—} & \(\rho_{e\to v}\) & — \\ \cline{3-6}  & & & & \multicolumn{2}{c|}{—} & \(\rho_{e\to v}\) & — \\ \cline{3-6}  & & & \multicolumn{2}{c|}{—} & \(\phi_{v}\) & \(y_{i}=b_{i}^{2}\) & \(\rho_{e\to g}\) & — \\ \cline{3-6}  & & & \multicolumn{2}{c|}{—} & \(\phi_{g}\) & \(h_{1}=\overline{\psi}\) & \(\overline{\psi}=\sum_{i}y_{i}\) \\ \cline{3-6}  & & & & \multicolumn{2}{c|}{—} & \(\rho_{v\to g}\) & — \\ \cline{3-6}  & & & & \multicolumn{2}{c|}{—} & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 3: Power Method as a Graph Network - Iterative Layers

\begin{table}
\begin{tabular}{|l l l||l|l|} \hline \multicolumn{3}{|c||}{} & \multicolumn{2}{c||}{Data} & \multicolumn{2}{c|}{Functions} \\ \multicolumn{3}{|c||}{Fixed} & \multicolumn{1}{c||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{2}{c|}{Aggregation} \\ \hline \hline \multirow{3}{*}{Edge} & \(A_{ij}\) & \(c_{ij}\) & \(\phi_{e}\) & \(c_{ij}=A_{ij}b_{j}\) & \(\rho_{e\to v}\) & \(\overline{c}_{i}=\sum_{j}c_{ij}\) \\  & — & \(b_{i},y_{i}\) & \(\phi_{v}\) & \(b_{i}=\overline{c}_{i}\) & \(\rho_{e\to g}\) & — \\ \cline{3-6}  & & & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \cline{3-6}  & & & \multicolumn{2}{c|}{—} & \multicolumn{2}{c|}{—} \\ \cline{3-6}  & & & \multicolumn{2}{c|}{—} & \(\phi_{e}\) & — \\ \cline{3-6}  & & & \(\phi_{v}\) & \(y_{i}=b_{i}^{2}\) & \(\rho_{e\to g}\) & — \\ \cline{3-6}  & & & \multicolumn{2}{c|}{—} & \(\rho_{v\to g}\) & — \\ \cline{3-6}  & & & \multicolumn{2}{c|}{—} & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 3: Power Method as a Graph Network - Rayleigh Quotient
expensive.

While Jacobi/Chebshev propagate information slowly across the mesh, it can be shown that these methods often efficiently produce approximate solutions where the error is smooth. The multigrid algorithm leverages this smoothing property in conjunction with a hierarchy of meshes \(\mathcal{G}_{\ell}\) such as those depicted in Figure 2.

One simple example of a multigrid iteration employs \(k\) Jacobi iterations on the \(A_{1}x_{1}=b_{1}\) system to generate an approximate solution \(x_{1}^{(k)}\). Since the error is now smooth, it can be accurately approximated on a coarser mesh. Thus, the multigrid idea is to improve the solution using a coarser mesh. To do this, a residual is computed and then projected on to the \(\mathcal{G}_{2}\) mesh within the hiearchy and used as the right hand side to a second linear system \(A_{2}x_{2}=b_{2}\). Here, the matrix \(A_{2}\) is a coarser discrete version of the Poisson equation constructed using \(\mathcal{G}_{2}\). The \(A_{2}\) matrix is a less accurate approximation to Poisson's equation, but it is only used to correct the approximation obtained on the finer mesh. One can apply \(k\) Jacobi iterations to this \(A_{2}\) system to produce \(x_{2}^{(k)}\). Notice that these iterations are less expensive than applying Jacobi to the \(A_{1}\) system. We can repeat this same process (apply Jacobi to \(A_{\ell+1}x_{\ell+1}=b_{\ell+1}\) where \(b_{\ell+1}\) is a projected version of the residual for the \(A_{\ell}\) system) for all meshes in the hierarchy. In each case, the approximate solution on the \(\ell^{th}\) mesh is a correction to the solution on the previous mesh system. To complete the AMG iteration, we must add together the individual corrections. That is, the approximate solution at the end of a single AMG iteration is given recursively by \(\tilde{x}_{\ell}=x_{\ell}^{(k)}+P_{\ell}\tilde{x}_{\ell+1}\) where \(P_{\ell}\) is a rectangular matrix that interpolates (or prolongates) approximate solutions associated with the \(\ell+1\) mesh to the \(\ell\) mesh. Here \(1\leq\ell\leq N_{levels}-1\) and on the coarsest mesh we take \(\tilde{x}_{N_{levels}}=x_{N_{levels}}^{(k)}\). Often a Gaussian elimination solver is used on the coarsest mesh as its cost is negligible when the matrix \(A_{N_{levels}}\) is sufficiently small. Notice that in this case, one AMG iteration propagates information from the clamps throughout the entire mesh. While there is some cost to the coarse level computations, this cost is small relative to the \(A_{1}\) computations if the coarser meshes are significantly smaller than the finest mesh. While it is generally difficult for application developers to generate the coarse operators \(\mathcal{G}_{\ell},A_{\ell}\), and \(P_{\ell}\) (\(P_{\ell}^{T}\) is commonly used for residual projections), an AMG algorithm automates this entire process. That is, all of the additional operators/meshes are generated within the AMG method using graph algorithms that coarsen the matrix graph, define _algebraic_ interpolation, and define the coarse \(A_{\ell}\)'s.

#### 3.6.1 Strength of Connection

Automatically coarsening a matrix graph is difficult. In our example, we need a coarse graph that can approximate errors that remain after \(k\) Jacobi steps. Several different AMG algorithms are used in practice each with its own coarsening procedure. Typically, a strength-of-connection algorithm is applied as a pre-cursor to coarsening. The basic idea is that some graph edges (or off-diagonal matrix nonzeros) may contribute very little to the information flow. For example, suppose one corner region of our plate consists of a near insulating material (e.g., extending from the corner to the closest segment of the circle). Edges between the conducting and insulating regions should be classified as weak so that they can effectively be ignored during the graph coarsening process. Unfortunately, improper edge classification can ruin AMG's impressive convergence properties and remains an active research topic that could possibly be improved by new ML algorithms. We present two common strength-of-connection algorithms that are based on evaluating the relative size of off-diagonal matrix entries. However, we note that these algorithms may not be appropriate for complex matrix systems. Unlike the previously discussed methods that output information at graph vertices, strength of connection algorithms output data on graph edges. Specifically, a strength of connection

Figure 2: Left: sample mesh hierarchy (clamps shown on \(\mathcal{G}_{1}\)). Right: data movement within an AMG iteration (Jacobi iterations occur at tan circles and Gaussian elimination occurs at dark brown circle).
algorithm produces a matrix \(S(A)\) which has the same sparsity pattern as the input matrix \(A\).

For the (non-symmetric) smoothed aggregation AMG strength of connection, this computation is,

\[S_{ij}=\frac{A_{ij}^{2}}{A_{ii}A_{jj}}. \tag{2}\]

The GNN layer version of this algorithm, which is shown in Table 8, relies on fixed matrix edge data, \(A_{ij}\), as well as fixed diagonal values, \(A_{ii}\), stored on the vertices. Note that this GNN allows for self edges, so the diagonal values are stored on both the vertices and the edges. The only calculation takes place in the edge update function, \(\phi_{e}\), which simply implements (2).

Similar to the smoothed aggregation strength of connection, the classical strength of connection algorithm provides output data on graph edges to construct a strength of connection matrix \(S(A)\) that has the same sparsity as the input matrix \(A\). For the classic strength of connection, this computation is,

\[S_{ij}=\frac{-A_{ij}}{\max_{k\neq i}\{-A_{ik}\}}. \tag{3}\]

The GNN representation of this algorithm, which is shown in Table 9, relies only on fixed matrix edge data, \(A_{ij}\). The maximum negative value on each row is first determined in the aggregation phase, \(\rho_{e\to v}\), and is then propagated to the vertices. Finally, the edge update is used to finish the calculation in (3).

One final step is needed to finish the classic strength of connection, and that is the dropping of weak connections. For a suitable \(0<\tau\leq 1\), drop the nonzero values where \(S_{ij}<\tau\). The sparsity pattern of \(S_{ij}\) then provides the classic strength of connection. This could be done as a post-processing step, or, to keep with the theme of utilizing ideas from graph neural networks, we can modify the final edge update to utilize a non-linear activation function, such as the **step** function:

\[\text{step}(x)=\begin{cases}1&x>0\\ 0&x\leq 0\end{cases}\]

to zero out the weak connections based on \(\tau\). This results in

\[\hat{S}_{ij}=\text{step}\left(\frac{-A_{ij}}{\max_{k\neq i}\{-A_{ik}\}}-\tau \right), \tag{4}\]

where the strong connections are the non-zero entries of \(\hat{S}\).

\begin{table}
\begin{tabular}{|l l l||l l|l l|} \hline \multicolumn{3}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\ \multicolumn{3}{|c||}{Fixed} & \multicolumn{1}{c||}{Mutable} & \multicolumn{3}{c|}{Updates} & \multicolumn{3}{c|}{Aggregation} \\ \hline \hline \multirow{3}{*}{Edge} & \multirow{3}{*}{\(A_{ij}\)} & \multirow{3}{*}{\(S_{ij}\) [Output]} & \multicolumn{3}{c|}{Layer 1} \\ \cline{3-3} \cline{5-8}  & & & \(\phi_{e}\) & — & \(\rho_{e\to v}\) & \(\overline{\epsilon}_{i}=\max_{j\neq i}\{-A_{ij}\}\) \\ \cline{3-3} \cline{5-8} Vertex & — & \(v_{i}\) & \(\phi_{v}\) & \(v_{i}=\overline{\epsilon}_{i}\) & \(\rho_{e\to g}\) & — \\ \cline{3-3} \cline{5-8} Global & — & — & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 9: Classic Strength of Connection as a Graph Network

\begin{table}
\begin{tabular}{|l l l||l l|l l|} \hline \multicolumn{3}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\ \multicolumn{3}{|c||}{Fixed} & \multicolumn{1}{c||}{Mutable} & \multicolumn{3}{c|}{Updates} & \multicolumn{3}{c|}{Aggregation} \\ \hline \hline \multirow{3}{*}{Edge} & \multirow{3}{*}{\(A_{ij}\)} & \multirow{3}{*}{\(S_{ij}\) [Output]} & \multicolumn{3}{c|}{Layer 1} \\ \cline{3-3} \cline{5-8}  & & & \(\phi_{e}\) & \(S_{ij}=A_{ij}^{2}/(A_{ii}A_{jj})\) & \(\rho_{e\to v}\) & — \\ \cline{3-3} \cline{5-8} Vertex & \(A_{ii}\) & — & \(\phi_{v}\) & — & \(\rho_{e\to g}\) & — \\ \cline{3-3} Global & — & — & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 8: Smoothed Aggregation Strength of Connection as a Graph Network
#### 3.6.2 Direct Interpolation

We present one of the simplest AMG interpolation schemes, though more advanced/robust algorithms are generally preferred. The interpolation algorithm provides output on the graph edges, similar to the strength of connection algorithms. Specifically, we seek the nonzeros entries of \(P_{\ell}\). In the description that follows, we omit the sub-script \(\ell\) to simplify notation. The algorithm relies on first partitioning the vertex set \(\mathcal{V}\) into \(\mathcal{F}\) and \(\mathcal{C}\) such that \(\mathcal{F}\cap\mathcal{C}=\emptyset\) and \(\mathcal{F}\cup\mathcal{C}=\mathcal{V}\) (cf., [7]). The F-vertices only exist on the fine grid while C-vertices exist on both fine and coarse grids. We assume that the rows/columns of A (and hence the vertices) are ordered so that all C-vertices are numbered before F-vertices. In addition to partitioning, we require a strength of connection matrix, \(S\), as discussed in Section 3.6.1.

The direct interpolation operator is derived from input matrix \(A\) such that if \(i\) is an F-vertex, then

\[P_{ij}=-A_{ij}\ \frac{\sum_{k\in N_{i}}A_{ik}}{A_{ii}\sum_{k\in C_{i}^{ \prime}}A_{ik}} \tag{13}\]

where \(N_{i}\) denotes the neighbors of vertex \(i\), and \(C_{i}^{\ast}\) denotes the strong, coarse neighbors of vertex \(i\). Specifically \(j\in C_{i}^{s}\) if and only if \(j\in N_{i}\), \(j\in\mathcal{C}\) and \(S_{ij}-\tau>0\). If \(i\) is instead a C-vertex,

\[P_{ij}=\begin{cases}1&j=i\\ 0&j\neq i\end{cases}.\]

A derivation of (13) can be found in [9]. One can see that the \(i^{th}\) row of \(P\) is simply a weighted sum where weights are proportional to \(A_{ij}\) and normalized by the fraction term, This fraction does not depend on \(j\) and guarantees that the sum of the nonzeros in the \(i^{th}\) row of \(P\) ( \(=\sum_{k\in C_{i}^{\prime}}P_{k}\) ) is equal to one when the sum of \(A\)'s nonzeros in the \(i^{th}\) row is zero. This is generally true when \(A\) represents differentiation as the derivative of a constant function is zero, which is effectively approximated by multiplying \(A\) by a vector where all entries are constant.

The GNN representation of this algorithm, which is shown in Table 10, requires two layers. We use the coarse/fine splitting information and the matrix diagonal as fixed vertex features. We also need a mutable, temporary vertex feature \(\alpha\) for communicating information between the two GN layers. As edge features, we use the off-diagonal entries of \(A\) and the strength of connection information as fixed features, and the output weight, \(P_{ij}\) is edge-specific and mutable.

As the algorithm proceeds, the edge update in the first layer is used to pass the coarse/fine splitting information,

\[C_{j}=\begin{cases}1&j\in\mathcal{C}\\ 0&j\in\mathcal{F}\end{cases},\]

of the receiving vertex to the edge. Next, the edge aggregation is used to perform the necessary summations since all the quantities are specific to the \(i\)th vertex. In the vertex update of the first layer, the \(A_{ii}\) factor is multiplied through. Finally, in the second layer, only the edge update function is utilized, which incorporates the \(A_{ij}\) value into the weight and also assigns zeroes to all rows associated with coarse points. After the computation of the GNN is complete, some additional post-processing is necessary. As mentioned previously, the GNN generates an updated graph with the same structure as the input graph (which is square since it is derived from \(A\) and \(S\)), but the prolongation operator is rectangular. To get a complete prolongation operator, two additional operations would be necessary: (a) setting the diagonal to one and (b) removing all columns associated with fine vertices. The first of these operations can be wrapped into the edge update in the second layer of the GNN network, but the second must be completed after the computation of the entire graph neural network. In doing so, the standard direct interpolation operator is obtained.

## 4 Training Within GNN Frameworks

The GNN models considered in the previous section used prescribed aggregation and update functions. Typically, update functions used are parameterized multi-layer perceptron networks, while the variadic aggregation functions are prescribed. This section gives two examples that illustrate how learned parameters can be introduced into update functions so that the GNN model can be trained to perform a numerical task.

Code AvailabilityThe code for the examples is this chapter is provided in the repository at [https://github.com/sandialabs/gnn-applied-linear-algebra/](https://github.com/sandialabs/gnn-applied-linear-algebra/). Only PyTorch code is provided. The repository includes code to create all datasets, implement the given GNNs, and train the models.

### Learning a Diagonal for the Jacobi Iteration

As discussed, an iteration of the weighted Jacobi method for solving \(Ax=b\) is defined by

\[x^{k+1}=x^{k}+\omega D^{-1}(b-Ax^{k}), \tag{14}\]
where the solution at the \(k^{th}\) Jacobi iteration is updated based on the associated residual. The \(i^{th}\) entry of the residual is scaled by \(\omega/D_{ii}\) to update the \(i^{th}\) solution entry. Here, \(\omega\) is a user-provided scaling value and \(D\) is a diagonal matrix whose nonzeros are given by \(D_{ii}=A_{ii}\). Ideally, \(\omega\) is chosen so that convergence is attained in the fewest number of iterations. While this ideal \(\omega\) can be determined by computing eigenvalues of the matrix \(D^{-1}A\) when \(A\) is a symmetric positive definite matrix, in practice it is often chosen in an ad hoc fashion.

We pursue an alternative machine learning approach that selects a diagonal relaxation operator \(\bar{D}\) based on the matrix. A generalized Jacobi method (or scaled Richard iteration) is given by

\[x^{k+1}=x^{k}+\bar{D}^{-1}(b-Ax^{k}), \tag{4.2}\]

where the diagonal matrix \(\bar{D}\) is not defined by \(A\)'s diagonal entries but is instead chosen by a machine learning algorithm to reduce the number of iterations required to reach convergence. We follow the framework introduced in [31] for learning the diagonal of the generalized Jacobi iterative method, though we consider a different target class of matrices to demonstrate the approach focusing on Jacobi as a relaxation method in multigrid. Within multigrid, the relaxation method should reduce high frequency errors as a complement to the coarse grid correction that addresses low frequency errors.

As an example, we discretize a 2D Poisson operator

\[\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\]

on the unit square \([0,1]\times[0,1]\) with homogeneous Dirichlet boundary conditions: \(u(x,0)=u(0,y)=u(x,1)=u(1,y)=0\). The domain is tiled using linear finite elements with uniformly shaped quadrilaterals, except for a small, thin band of tall and skinny quadrilaterals, such as the mesh shown in Figure 4.1. Note that the horizontal location of the band and its width can be changed to produce different variations of the problem. This example is meant to highlight a situation where a generalized \(\bar{D}\) matrix within Jacobi's method might be advantageous, and should not be misconstrued as good meshing practice.

In the matrix arising from discretizing this system, the abrupt change in mesh spacing yields significant differences when comparing rows associated with mesh nodes far from the band versus mesh nodes within the band. Two matrix stencils that highlight this variation are

\[\frac{1}{6}\left[\begin{array}{ccc}-2&-2&-2\\ -2&16&-2\\ -2&-2&-2\end{array}\right]\qquad\quad\text{and}\qquad\quad\quad\frac{\beta}{6 h}\left[\begin{array}{ccc}-1-\left(\frac{h}{\beta}\right)^{2}&-4+2\left(\frac{h}{ \beta}\right)^{2}&-1-\left(\frac{h}{\beta}\right)^{2}\\ 2-4\left(\frac{h}{\beta}\right)^{2}&8+8\left(\frac{h}{\beta}\right)^{2}&2-4 \left(\frac{h}{\beta}\right)^{2}\\ -1-\left(\frac{h}{\beta}\right)^{2}&-4+2\left(\frac{h}{\beta}\right)^{2}&-1- \left(\frac{h}{\beta}\right)^{2}\end{array}\right].\]

The left stencil corresponds to a mesh node like the node surrounded by a square in Figure 4.1 which is away from the boundary and the thin band (so locally the mesh is uniformly spaced in both the \(x\) and \(y\) directions). The right stencil corresponds to a mesh node like the node surrounded by a circle in Figure 4.1 which in the thin band and away from the boundary; the vertical spacing is \(h\) and the horizontal spacing (both to the left and right of the node) is \(\beta\). When \(h/\beta=1\), the left and right stencils correspond. Notice that for \(\beta\ll h\), the \(h/\beta\) terms dominate and one can see that the left and right stencils are very different. While our example is artificially devised to stress Jacobi's method, this type of

\begin{table}
\begin{tabular}{|l l l l||l l l|} \hline \multicolumn{4}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\  & Fixed & Mutable & Updates & Aggregation \\ \hline \hline \multirow{4}{*}{Edge} & \multirow{4}{*}{\(A_{ij},\hat{S}_{ij}\)} & \multirow{4}{*}{\(v_{ij},P_{ij}\) [Output]} & \multicolumn{4}{c|}{Layer 1} \\ \cline{3-6}  & & & & \(\phi_{e}\) & \(v_{ij}=C_{j}\) & \(\rho_{e\to v}\) & \(\overline{\gamma}_{i}=\frac{\sum_{j}A_{ij}}{\sum_{j}A_{ij}v_{ij}\hat{S}_{ij}}\) \\ \cline{3-6}  & & & & \(\phi_{v}\) & \(\alpha_{i}=\frac{1}{A_{ii}}\overline{\gamma}_{i}\) & \(\rho_{e\to g}\) & — \\ \cline{3-6}  & & & & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \cline{3-6}  & & & & \multicolumn{3}{c|}{Layer 2} \\ \cline{3-6}  & & & & \(\phi_{e}\) & \(P_{ij}=(1-C_{i})(-A_{ij}\alpha_{i})\) & \(\rho_{e\to v}\) & — \\ \cline{3-6}  & & & \(\phi_{v}\) & — & \(\rho_{e\to g}\) & — \\ \cline{3-6}  & & & & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 10: Direct Interpolation Kernel as a Graph Network
abrupt stencil change occurs in more realistic scenarios. For instance in [10], generalized Jacobi iterations are used for a discontinuous Galerkin discretization where penalties glue meshes together.

#### 4.1.1 Datasets

The training, test and validation datasets are generated by varying the location and width of the band in the mesh. The steps for producing the datasets are as follows:

1. Create a 2D uniform mesh which contains \(N_{y}\) points in each direction
2. Randomly select \(\beta\) between \(\beta_{\min}\) and \(\beta_{\max}\)
3. Randomly select an existing \(x\)-coord from the mesh in step 1
4. Build the band mesh by placing new mesh vertices a distance of \(\beta\) to the left and right of all vertices with the \(x\)-coordinate in step 3.
5. Build the matrix based on a finite element discretization
6. Repeat steps 1-5 to generate 1000 matrices
7. Matrices 1-800 become the training dataset
8. Matrices 801-850 become the validation set
9. Matrices 851-1000 become the test set

The training set defines the loss function used by gradient descent to determine model parameters that yield a reduced loss value. Using the validation set, the model (built with the training set) can be evaluated over a range of different model and algorithm choices (referred to as hyperparameters) such as model architectures and optimizers. The test dataset is used to evaluate the final model obtained after all hyperparameter tuning has occurred. The test set is never used during the model selection process.

#### 4.1.2 Input/Output Attributes

The input edge attributes are the nonzero matrix entries \(A_{ij}\) excluding the diagonal. The matrix diagonal, \(A_{ii}\) is included as a a vertex attribute for each \(i\). The output vertex attribute will be \(d_{i}=\bar{D}_{ii}^{-1}\) for the generalized Jacobi method. See the _Data_ section of Table 1.

#### 4.1.3 Loss Function

A neural network is trained by identifying parameter values that approximately minimize or at least significantly reduce the value of a loss (also known as an objective) function. One standard machine learning technique involves finding _optimal_ answers (typically from an expensive process) and then training the neural network to match these optimal values; a technique referred to as supervised learning. In our case, this would require computing a set of optimal \(\bar{D}_{ii}\) values for each training scenario using a (possibly expensive) numerical procedure. Then, a loss function would be defined that minimizes the difference between these pre-computed values and the GNN output version of

Figure 1: An example mesh with a thin band of elements used for training and testing. Here \(h=\frac{1}{l}\) and \(\beta=0.05\)
these values. Unfortunately, a tractable numerical procedure for pre-computing the \(\bar{D}_{ii}\) values is not apparent given the fact that numerous eigenvalue calculations would be needed to determine an _optimal_ high dimensional vector (in \(\mathbb{R}^{N}\)) to define \(\bar{D}\) for each training case.

Instead we define the loss function based on the desired performance of the Jacobi method itself. We would like to minimize the number of required Jacobi iterations to reach some specified convergence tolerance. However, efficient optimization methods for neural network training require computing the gradient of the loss function and since the number of iterations is a discrete quantity, we cannot calculate a gradient with respect to it. Thus, the number of iterations cannot be used and a differentiable quantity that measures the performance of a Jacobi iteration will be used instead. From the theory of iterative methods, we that the asymptotic convergence rate of a linear iterative method is given by the spectral radius of the error propagation matrix. The error propagation matrix for generalized Jacobi is \(I\!-\!\bar{D}^{-1}A\), which is obtained by substituting \(e^{k}=x^{k}\!-\!A^{-1}b\) and \(e^{k+1}=x^{k+1}\!-\!A^{-1}b\) in (4.2). Thus, we seek to minimize the spectral radius of this \(I\!-\!\bar{D}^{-1}A\) matrix.

Here, we face another obstacle: standard automatic differentiation tools do not support eigenvalue algorithms. Even if support was added, [34] shows that standard eigenvalue decomposition algorithms and power method algorithms give numerically unstable gradients. Thus, we need a different method to approximate the spectral radius of a matrix which gives a numerically stable eigenvalue estimate. In [31], it is proved that the spectral radius of a matrix, \(B\) can be approximated by

\[\max_{i=1,2,...,m}\{\|B^{K}\hat{u}_{i}\|^{\frac{1}{K}}\}\]

where \(K\) is a user defined number of iterations, \(\{\hat{u}_{i}\}_{i=1}^{m}\) is a set of vectors randomly chosen from the surface of the unit sphere. For multigrid relaxation, we instead want to maximize the performance of Jacobi with regard to only high-frequency errors. We can characterize the high frequency space using a discrete sine transformation matrix \(V\) that is defined such that each column has the general form

\[\alpha\sin(\theta_{x}\pi x)\sin(\theta_{y}\pi y).\]

Specifically, each of the \(N_{y}^{2}\) columns corresponds to a unique \((\theta_{x},\theta_{y})\) pair chosen from \(\theta_{x}=1,...,N_{y}\) and \(\theta_{y}=1,...,N_{y}\). The scalar \(\alpha\) is picked to ensure that the norm of each column is one and for the results in Section 4.1.6, we use \(N_{y}=38\). Partitioning \(V\) into low frequency and high frequency columns, we have \(V=[V_{lf}\;\;V_{hf}]\) where the low frequency modes correspond to pairs where both \(\theta_{x}\leq N_{y}/2\) and \(\theta_{y}\leq N_{y}/2\) while all other modes define \(V_{hf}\). As we seek to minimize high frequency errors via Jacobi relaxation, for each matrix \(A^{(j)}\) we train the GNN to find \(\bar{D}\) that minimizes

\[\max_{i=1,2,...,m}\{\|(I-\bar{D}^{-1}A^{(j)})^{K}\hat{u}_{i}\|^{\frac{1}{K}}\} \tag{4.3}\]

where each \(\hat{u}_{i}\) is a randomly chosen column of \(V_{hf}\). Therefore, if \(N\) is the number of matrices in the dataset being evaluated (training, validation, or test), the overall loss function is

\[\mathcal{L}=\sum_{j=1}^{N}\max_{i=1,2,...,m}\{\|(I-\bar{D}^{-1}A^{(j)})^{K}\hat {u}_{i}\|^{\frac{1}{K}}\} \tag{4.4}\]

For the results given in Section 4.1.6, we use \(K=3\) and \(m=20\).

#### 4.1.4 Architecture

In general the choice of architecture requires a search over range of parameters known as an ablation study. However, for this example, a fixed architecture for the neural network is chosen as a simple demonstration. The GNN consists of a single GNN layer with a vertex update, as can be seen in Table 1. The edge to vertex aggregation function, \(\rho_{e\to v}\), calculates the min, mean, sum and max of the edge attributes, resulting in four quantities per vertex. These four attributes are concatenated with the single vertex attribute (\(A_{ii}\)) to form a vector in \(\mathbb{R}^{5}\) input into the vertex update function. The vertex update function, \(\phi_{v}\) is a feed-forward neural network consisting of 3 layers. For more information about the architecture of this GNN, see Appendix A.1.

\begin{table}
\begin{tabular}{|l l l||l|l|l|} \hline \multicolumn{3}{|c||}{Data} & \multicolumn{3}{c|}{Functions} \\  & Fixed & Mutable & \multicolumn{1}{c|}{Updates} & \multicolumn{1}{c|}{Aggregation} \\ \hline \hline Edge & \(A_{ij}\) & — & \multicolumn{3}{c|}{Layer 1} \\ \cline{4-6} Vertex & \(A_{ii}\) & \(d_{i}\) [Output] & \(\phi_{e}\) & — & \(\rho_{e\to v}\) & \(\overline{c}=[min,mean,sum,max]\) \\ Global & — & — & \(\phi_{g}\) & — & \(\rho_{v\to g}\) & — \\ \hline \end{tabular}
\end{table}
Table 1: Training a Diagonal for Jacobi as a Graph Network
#### 4.1.5 Training Methodology

The model is implemented in PyTorch [25], using the PyTorch Geometric library [12]. The GNN is trained for 62 epochs using the Adam optimization method [21] with batch size of 100. The number of epochs can have a large influence on the overall success of the training phase. Too few epochs and the training may not reach peak performance, while too many epochs could lead to the model "over-fitting". In this case, the model fits the data too tightly and so the ability of the model to generalize to new data degrades. There are different methods to address over-fitting; we choose to stop training when the performance of the model on the validation set starts to degrade. For this experiment, the number of epochs is selected by training the model for 100 epochs, but stopping the training when the loss on the validation set is minimized. By stopping when the validation error is minimized, we obtain good generalization properties for future data.

#### 4.1.6 Results

We compare the result of the generalized Jacobi iteration using the learned model \(\bar{D}\) with a traditional weighted Jacobi using three different standard weights: standard Jacobi (\(\omega=1\)), a "default" weight (\(\omega=\frac{2}{3}\)), and a "classical optimal" weight (\(\omega_{\text{co}}=\frac{2}{\lambda_{\min}+\lambda_{\max}}\) where \(\lambda_{\min},\lambda_{\max}\) are the minimum and maximum value of \(D^{-1}A\) respectively). Note that the "classical optimal" weight is only optimal when Jacobi is used directly to solve \(Ax=b\) as opposed to Jacobi relaxation within a multigrid algorithm. To compare the different Jacobi relaxation procedures, we compute the maximum eigenvalue(s) of \(I-\omega V_{th}^{T}D^{-1}AV_{hf}\) (or of \(I-V_{th}^{T}\bar{D}^{-1}AV_{hf}\) for the generalized Jacobi method), thereby investigating each method's ability to damp high frequency error. The 10 largest eigenvalues associated with two of the test matrices are shown in Figure 2.

We can see from the two sample matrices that in some cases the learned diagonal outperforms the rest of the methods, while in some cases, the \(\omega_{\text{co}}\) weight performs the best. It would therefore be helpful to see a full comparison across all test matrices. In Figure 3, we compare the distributions of maximal eigenvalues of the standard methods with respect to the learned diagonal method. Hence, in the figure, all counts left of the 0.0 line indicate matrices where the associated method performed better than the learned method, while all counts to the right of the

Figure 3: Distributions of the differences in the maximum eigenvalues from the learned diagonal and the given constant weights for the Jacobi method. Counts left of the 0.0 line indicate matrices for which the associated method has lower maximum eigenvalue than the learned method (and by how much), while counts right of the 0.0 line indicate matrices for which the associated method has larger maximum eigenvalue than the learned method.

Figure 2: Top 10 largest eigenvalues of the generalized Jacobi relaxation method for two different test matrices, comparing the different methods of selecting the \(\omega\) parameter in weighted Jacobi to the learned diagonal from the GNN.
\(0.0\) line indicate matrices where the associated method performed worse than the learned method. From these results, we see that the learned method gives a better result than the standard non-weighted and \(\omega=\frac{2}{3}\) weighted Jacobi method in all cases and outperforms the \(\omega_{co}\) weight in more than \(75\%\) of the test matrices.

The shape of the right-most histogram in Figure 3 suggests that there might be areas of the test space where the learned method performs better and areas where the \(\omega_{co}\) performs better. To discover what patterns might exist, we map test matrices with regard to the width and location of the band in Figure 4 and color their associated dot with the top performing method.

We can see from this plot that there is a clear pattern with regard to which matrices are more effectively solved using the learned diagonal versus the \(\omega_{co}\), constant weight. This makes intuitive sense because for larger widths, the mesh is "more uniform" and thus the correctly chosen constant will perform well, while for smaller widths, allowing the weight to vary on a per-row basis allows for more flexibility to address the more drastically different elements.

Finally, comparison plots for the scaled diagonal values arising from the \(\omega_{co}\) method and the learned diagonal are given in Figure 5 for a test matrix with \(h=0.003\) and where the band is approximately located at \(x=0.405\).

There are several observations which can be made from the images. First, we observe that in general, the learned method yields larger diagonal values than the \(\omega_{co}\) method. It can also be observed that the learned method produces a larger diagonal in the center of the band than the \(\omega_{co}\) method. Most drastically, we observe in the learned method that boundary points (especially the corners) have larger diagonal values than the interior points. This is true even in the band, where the boundary points are also larger than the interior points. This is a modification that the traditional weighted Jacobi method cannot replicate since the diagonal of \(A\) is the same at the boundaries and interior.

Figure 4: A plot of all the test matrices according to the width and location of their band, where color represents the winning method. Red represents the learned diagonal method while green represents the \(\omega_{co}\) constant weight method.

Figure 5: Diagonals at each mesh point for the \(\omega_{co}\) method and the learned method
The shapes in Figure 5 demonstrate why the learned method is able to outperform the "classical optimal" weighting for Jacobi's method. For the standard weighted Jacobi, only a constant multiple of \(D^{-1}A\) is allowed, which forces the shape of the diagonal to match the shape of the inverse of the diagonal of \(A\). This shape cannot be adjusted, nor can the relative distances be changed (only scaled). In our new learned method, no such limitation exists, hence the space from which the learned method can select the diagonal values is much broader. This broader space presents greater opportunity to find an optimal diagonal.

### Example: Training a GNN to Determine Diffusion Coefficients

The strength of connection algorithms discussed in Section 3.6.1 are heuristic schemes with several shortcomings. One could instead consider GNN models that learn a more sophisticated measure of strength. In this section, we consider the related task of determining diffusion coefficients given matrix stencil and coordinate information. Though not equivalent to determining strength measures, there is a rough connection between the relative sizes of diffusion coefficients in different directions and the strength of the matrix connections aligned with these directions.

Consider the diffusion operator

\[\sum_{i=1}^{2}\sum_{j=1}^{2}\frac{\partial}{\partial x_{i}}\left[D_{ij}(x,y) \frac{\partial}{\partial x_{j}}\right] \tag{10}\]

defined on the unit square with periodic boundary conditions. \(D(u,x,y)\) is assumed to have the form

\[D(x,y)=\begin{bmatrix}\alpha(x,y)&0\\ 0&\beta(x,y)\end{bmatrix}.\]

The goal of this learning task is to predict the diffusion coefficients \(\alpha(x,y)\) and \(\beta(x,y)\) at each grid point using information from the coefficient matrix, coordinates of vertices, and mesh spacing.

#### 4.2.1 Datasets

All the datasets consist of data generated from discretizing (10) with the finite element method. Using the following six-step procedure, 1000 matrices are generated:

1. Select a random integer \(N\in[80,100]\) and let the mesh resolution be \(h=\frac{1}{N}\)
2. Select \(\theta_{\alpha,x},\theta_{\alpha,y},\theta_{\beta,x},\theta_{\beta,y}\) each with uniform probability from \(\{i:i\in\mathbb{Z}\text{ and }0\leq i\leq 6\}\)
3. Define \(\alpha(x,y)=\cos\left(\theta_{\alpha,x}\pi x\right)^{2}\cos\left(\theta_{ \alpha,y}\pi y\right)^{2}\)
4. Define \(\beta(x,y)=\cos\left(\theta_{\beta,x}\pi x\right)^{2}\cos\left(\theta_{\beta, y}\pi y\right)^{2}\)
5. Discretize (10) on a uniform 2D quadrilateral mesh with resolution \(h\) in both \(x\) and \(y\) directions to construct the matrix operator \(A\)
6. Generate input features and output targets for this matrix

Notice that diffusion coefficients are chosen so that there is no discontinuity over the periodic boundary conditions. Thus, the diffusion fields vary smoothly over the domain.

The data set generated is divided into a training set (700 matrices), validation set (200 matrices), and test set (100 matrices). The training set is used to determine model parameters that yield a sufficiently small loss. The validation set is used to determine the model architecture, as will be explained shortly. Finally, the test set is only used once at the end of the study to evaluate the final model's effectiveness on unseen data. The test set is evaluated in the results sections where it is used as a trustworthy indication of how the chosen model will perform on new, completely unseen data.

#### 4.2.2 Input/Output Attributes

The left column of Table 2 describes the input/output attributes. The vertex attributes \(v_{i}\) consist of the matrix diagonal entries. That is, \(v_{i}=A_{ii}\). The edge attributes, \(c_{ij}\) are constructed as \(c_{ij}=(A_{ij},x_{\text{rel}},y_{\text{rel}})\) where \(x_{\text{rel}},y_{\text{rel}}\) are the relative differences between vertex \(i\) and vertex \(j\) in the \(x\) and \(y\) directions respectively, scaled by \(\frac{1}{h}\). For example, if vertex \(j\) is the vertex southeast of vertex \(i\), then the edge features for the edge from vertex \(i\) to vertex \(j\) are \(c_{ij}=(A_{ij},1,-1)\). There is a global attribute in this example, which is the mesh resolution: \(g=h\). The output of the neural network is the updated vertex attributes which are the predicted values for \(\alpha\) and \(\beta\) at each mesh vertex.

#### 4.2.3 Loss Function

The loss function is given by the mean squared error (MSE)

\[\mathcal{L}_{\omega}^{(k)}=\frac{1}{2N^{2}}\sum_{i=1}^{N^{2}}||\alpha^{(k)}(x _{i},y_{i})-\tilde{\alpha}_{\omega}^{(k)}(x_{i},y_{i})||_{2}^{2}+||\beta^{(k)} (x_{i},y_{i})-\tilde{\beta}_{\omega}^{(k)}(x_{i},y_{i})||_{2}^{2}\]

where the \((k)\) superscript denotes the matrix from the training set; \(\omega\) refers to the model's trainable parameters (determined by numerical optimization during training); the \((x_{i},y_{i})\) correspond to different mesh points and \(\tilde{\alpha}_{\omega}^{(k)}\) and \(\tilde{\beta}_{\omega}^{(k)}\) represents the GNN model predictions using the model parameter \(\omega\).
#### 4.2.4 Architecture

The model architecture employs an encoder, followed by the GNN outlined in Table 2. The encoder is applied to the edge and vertex attributes separately before the execution of the graph neural network. Specifically, an MLP, \(\text{Encoder}_{e}\), is applied to each set of edge attributes. Similarly, an MLP, \(\text{Encoder}_{v}\), is applied to each set of vertex attributes. Finally, a third MLP, \(\text{Encoder}_{g}\) is applied to the global attributes. Details on the design of the encoder MLPs can be found in Appendix A.2. Generically, encoders and decoders are common in GNNs. These are used to enhance the expressiveness of user input attributes with the aim to improve predictions and network trainability. Within the graph neural network layer, MLP neural networks are used for the edge update and vertex update functions.

This architecture was chosen after performing experiments using the validation set. The validation set allows researchers to experiment with different model parameters, architectures, optimizers, etc. without exposing the model to the test set prematurely. In our case, we investigated the following variations

* Number of GNN layers: 1, 2, or 3
* Number of MLP layers in the update functions: 1, 2, 3, or 4
* Width of the hidden layers in the MLP update functions: 8, 16, 32, or 64
* Encoder architecture: No encoder, 1 layer with 16 neurons, or 3 layers with 16 neurons
* Decoder architecture: No decoder, 1 layer with 16 neurons, or 3 layers with 16 neurons.

This yields a total of \(3\!\cdot\!4\!\cdot\!4\!\cdot\!3\!\cdot\!3=432\) total architectures combinations. To improve computational efficiency, we evaluated the different models in two stages. In the first stage, we employ \(25\%\) of the training set and \(25\%\) of the validation set to quickly identify the highest performing models among the \(432\) possibilities. The five best-performing models are then evaluated in stage two on the full validation set after they have been optimized using the full training set. As described in more detail in Appendix A.2, the best performing model of the \(432\) is parameterized by \(14,002\) trainable parameters (1 GNN layer, 2 MLP layers of width 16 in the update function and 3 layer 16 neuron encoders and decoders).

#### 4.2.5 Training Methodology

The model is again implemented in PyTorch [25], using the PyTorch Geometric library [12]. The Adam optimization method [21] is utilized to train the neural network weights with a batch size of 10. Of particular importance is the choice of the number of training epochs. As noted earlier, too few epochs may lead to poor accuracy while too many epochs may lead to overfitting where the model performs poorly on new, unseen data. We again choose to stop training when the performance of the model on the validation set starts to degrade. In our case, the minimum validation loss occurs after training the model for \(187\) epochs.

#### 4.2.6 Results

Figure 6 shows the MSE loss for the training and validation sets as a function of number of epochs. Each epoch is a forward and backward propagation to compute the gradient over all batches in the training set. Thus with a training set of \(700\) entries, and a batch size of \(10\) that is \(70\) gradient computations per epoch. The loss for both sets is attains a similar level (\(5.73\times 10^{-4}\) for validation and \(5.74\times 10^{-4}\) for testing). This similarity indicates that it is unlikely that over-training has occurred. Over-training is usually marked by a descending training error while simultaneously increasing validation error. This indicates the network does not generalize to data outside of the training set.

We also consider model performance in relation to the frequency of the diffusion functions. To study this behavior, assume \(D(x,y)\) takes the form

\[D(x,y)=\begin{bmatrix}\alpha(x,y)&0\\ 0&\alpha(x,y)\end{bmatrix}\]

where \(\alpha(x,y)\) has the form described in Section 4.2.1. Now, we can test the model on problems where \(\theta_{\alpha,x}\) and \(\theta_{\alpha,y}\) are selected in a grid from the set \(\{i:i\in\mathbb{Z}\text{ and }0\leq i\leq 16\}\). The frequency versus mean-squared error plot is given in Figure 7. Recall that the model is only trained for \(0\leq i\leq 6\), so

\begin{table}
\begin{tabular}{|l l l||l l|l l|} \hline \multicolumn{3}{|c||}{Data} & \multicolumn{4}{c|}{Functions} \\  & \multicolumn{1}{c}{Fixed} & \multicolumn{1}{c||}{Mutable} & \multicolumn{2}{c|}{Updates} & \multicolumn{2}{c|}{Aggregation} \\ \hline \hline \multirow{4}{*}{Edge} & — & \multirow{4}{*}{\(c_{ij}\) [Output]} & \multirow{4}{*}{\(c_{ij}\) [Output]} & \(\phi_{e}\) & \(c_{ij}=\text{Encoder}_{e}(c_{ij})\) & \(\rho_{e\to v}\) — \\  & & & \(\phi_{e}\) & \(v_{i}\) & \(\text{Encoder}_{v}(v_{i})\) & \(\rho_{e\to g}\) — \\ \cline{1-1}  & & & \(\phi_{g}\) & \(g\) & \(=\text{Encoder}_{g}(g)\) & \(\rho_{v\to g}\) — \\ \cline{1-1}  & & & \multicolumn{4}{c|}{Graph Neural Network} \\ \cline{5-6}  & & & \(\phi_{e}\) & \(c_{ij}=\textit{NN}_{e}(c_{ij},v_{i},v_{j},g)\) & \(\rho_{e\to v}\) & \(\overline{c}_{i}=[min,mean,sum,max]\) \\ \cline{1-1}  & & & \(\phi_{v}\) & \(v_{i}\) & \(=\textit{NN}_{v}(v_{i},\overline{c}_{i},g)\) & \(\rho_{e\to g}\) — \\ \cline{1-1}  & & & \(\phi_{g}\) & — & \(\rho_{v\to g}\) — \\ \hline \end{tabular}
\end{table}
Table 2: Training a Graph Network to Determine Diffusion Coefficients
all frequency combinations outside this interval are being extrapolated by the model. This portion of the subdomain is indicated by the shaded region. All matrices used in this study are generated using the same finite elements as previously, and all have 100 nodes in each direction yielding 10,000 by 10,000 matrix problems. The figure indicates that coefficients with low frequencies, those in the training set, are well approximated. Deviating from the training set, there are two sources of potential error, both of which contribute to the increase in the error observed in the figure. First, higher frequency coefficients suffer increasing error on a fixed (\(100\times 100\)) mesh as the number of points per wavelength decreases. The second source of error comes from the model extrapolating outside of the training set. In this context, examining the error away from the training frequencies we see that the error remains relatively small even with larger departures from the training data.

Connecting this result to strength-of-connection metrics, we remark that in general, it can be difficult to correctly classify strong and weak connections when \(\alpha\) is small and \(\beta\) is not small (or vice versa) when linear finite elements are used on quadrilateral meshes. Such a stencil where \(\alpha(x,y)=0.001\) and

Figure 4.6: Loss plots for training diffusion coefficients

Figure 4.7: Frequency vs mean-squared error for the trained model predicting point-wise diffusion coefficients from the finite element problem matrix and relative coordinates. The shaded region demonstrates the area where the training data was taken from, while the non-shaded region demonstrates where the model is extrapolating.
\(\beta(x,y)=0.8\) is shown below:

\[\left[\begin{array}{ccc}-0.1335&-0.533&-0.1335\\ \\ 0.266&1.068&0.266\\ \\ -0.1335&-0.533&-0.1335\end{array}\right]\]

Notice that if the following standard strength-of-connection metric is used:

\[|A_{ij}|\geq\theta\max_{k\neq i}|A_{i,k}|\Leftrightarrow A_{ij}\text{ is strong}\]

then \(\theta<0.25\) would (incorrectly) classify all four cardinal directions as strong. Allowing the trained model to predict the \(\alpha\) and \(\beta\) for this case, we obtain the following values:

\[\alpha=-0.00052847\quad\beta=0.71802.\]

As such, the model has correctly predicted this stencil as having small \(\alpha\) and with a reasonable approximation for \(\beta\). This provides an opportunity for giving better guidance to direct the multigrid coarsening procedure. Notice that the \(\alpha\) prediction is slightly negative. This is a by-product of using a LeakyReLU for the final activation allowing for negative outputs. The standard ReLU, while strictly positive, yields several training cases that have smaller derivative values, slowing down learning, and ultimately negatively impacting model accuracy. Applying a thresholding similar to a ReLU as a post-training step is possible, and would act as a more "reasonable" filter for outputs, but since we wished to evaluate the performance of our model, no additional steps were applied to the output of the GNN. This discrepancy also points to the possibility that more fine-tuning could improve the trained model, but as this is a more educational paper in nature, we forego adding any additional complexities to the model in this study.

## 5 Conclusion

We have examined graph neural networks (GNNs) from the perspective of numerical linear algebra. Our objective has been to highlight the close relationship between several sparse matrix tasks and GNNs. This relationship should not be surprising given the strong connections between sparse matrices and graphs. Several traditional numerical algorithms have been recast using a GNN representation. The GNN algorithms presented here will often be much less efficient than specialized sparse linear algebra libraries like Trilinos [18] or PETSc [3], but are instead intended to familiarize the reader with the inner workings and operations of GNNs by relating well-known algorithms to a GNN counterpart. We believe that GNNs may prove useful for many sophisticated linear algebra tasks where traditional algorithms have significant shortcomings. For example, one frequently recurring linear algebra theme is that of matrix approximation. The details vary by context, but the general idea is to find an inexpensive surrogate matrix that approximates the behavior of a large sparse matrix. For multigrid, this inexpensive matrix might be needed to approximate the action of the original matrix on a low frequency space. In other situations (e.g., reduced order modeling) an approximation of matrix vector products might be needed only for vectors that lie in a relatively small subspace driven by the simulation. These inexpensive approximations leverage smaller, sparser or lower rank matrices. Machine learning automates the task of generating approximation using data in cases that can be difficult with traditional approaches. The power of GNNs within machine learning is that their flexibility can address general sparse matrices. In this paper, we give two GNN learning examples to demonstrate, in simple cases, how GNNs can address some linear algebra tasks that can be incorporated into more traditional methods. While these examples are basic demonstrations of where GNNs can be used, we hope they will provide inspiration for how to incorporate GNNs into more complex problems.
