# Operation-Level Performance Benchmarking of Graph Neural Networks for
Scientific Applications

Ryien Hosseini1, 2, Filippo Simini1, Venkatram Vishwanath1,
{rhosseini, fsimini<EMAIL_ADDRESS>1 Argonne Leadership Computing Facility,
Argonne National Laboratory, Lemont IL USA  2 Department of Electrical
Engineering and Computer Science, University of Michigan, Ann Arbor MI USA

###### Abstract

As Graph Neural Networks (GNNs) increase in popularity for scientific machine
learning, their training and inference efficiency is becoming increasingly
critical. Additionally, the deep learning field as a whole is trending towards
wider and deeper networks, and ever increasing data sizes, to the point where
hard hardware bottlenecks are often encountered. Emerging specialty hardware
platforms provide an exciting solution to this problem. In this paper, we
systematically profile and select low-level operations pertinent to GNNs for
scientific computing implemented in the Pytorch Geometric software framework.
These are then rigorously benchmarked on NVIDIA A100 GPUs for several various
combinations of input values, including tensor sparsity. We then analyze these
results for each operation. At a high level, we conclude that on NVIDIA
systems: (1) confounding bottlenecks such as memory inefficiency often
dominate runtime costs moreso than data sparsity alone, (2) native Pytorch
operations are often as or more competitive than their Pytorch Geometric
equivalents, especially at low to moderate levels of input data sparsity, and
(3) many operations central to state-of-the-art GNN architectures have little
to no optimization for sparsity. We hope that these results serve as a
baseline for those developing these operations on specialized hardware and
that our subsequent analysis helps to facilitate future software and hardware
based optimizations of these operations and thus scalable GNN performance as a
whole.

##  I Introduction

In recent years, substantial research has been devoted to the development of
Graph Neural Networks (GNNs) [1]. GNNs have demonstrated to be powerful
machine learning tools, which are able to learn robust representations of data
in the form of graphs, and other graph-like data structures, such as point
clouds, natural networks, and manifolds [1, 2]. GNNs have applications in
several diverse scientific tasks. For example, molecules, proteins, social
networks, and knowledge networks can all be represented as graphs, and GNNs
have been shown to accurately and efficiently solve many tasks in these
applications [1], [2]. As the prominence of GNNs increase, so does the
importance of its runtime performance and ability to extend such networks to
larger data and model sizes.

In this paper, we profile and analyze state-of-the-art GNN models used in
several scientific applications, identify common operations of such models
when implemented in Pytorch/Pytorch Geometric [3], [4], and then systemically
benchmark these operations on a single NVIDIA A100 GPU. We evaluate operation
benchmarks on several dimensions, such as input sizes and data sparsity. Thus,
the primary purpose of this research is threefold:

  1. 1.

Aid software and hardware developers in identifying performance bottlenecks
and thus focus effort on future improvements in GNN development tools,

  2. 2.

Provide machine learning scientists performance insights so as to optimize
future GNN architectures

  3. 3.

Provide performance information on state-of-the-art hardware platforms for GNN
inference.

###  I-A Graph Neural Networks

At a high level, graph neural networks extend common neural network
operations, such as convolution, to non-euclidean data [1]. While traditional
deep learning applications such as computer vision operate on euclidean data,
the generalization of common deep learning operations enable to use of these
data-driven methods to solve a variety of scientific problems, in areas such
as particle simulation, protein folding, traffic prediction, social-network
related prediction tasks [1].

At the heart of GNNs is the message passing operation, wherein neighboring
nodes of a graph are iteratively updated based on differentiable and
permutation invariant functions [1] [4]. Recently, a plethora of powerful GNN
architectures have been developed. Indeed, many common deep learning
architectures, such as recurrent networks, attention based methods, and
variational autoencoders now have GNN equivalents.

A principal advantage of GNNs is their relative computational efficiencies.
Since the message passing algorithm operates on graphs, sparse inputs often
require far less computation than their equivalent traditional DNNs on their
euclidean-mapped counterparts. Thus, new operations that perform common
operations (e.g. matrix multiplication, transpose, etc.) on sparsely formatted
data have been introduced in many deep learning frameworks, such as Pytorch
[3] [4].

###  I-B Explicit Software Support for Sparse Operations

We focus this work on the Pytorch [3] deep learning framework. Recently, as
research interest in GNNs has increased, libraries such as Pytorch Geometric
[4] have been developed in order to provide explicit support for creating GNNs
and, by extension, sparse operations. Specifically, Pytorch Geometric is built
directly on Pytorch and provides support for a diverse set of GNN development
tasks, including end-to-end network implementations, efficient data handling
operations, and optimizations of sparse operations that are often used by
GNNs.

Pytorch Geometric’s operation level dependencies are torch_sparse,
torch_scatter, torch_cluster, and torch_spline_conv libraries. These libraries
provide optimizations over existing Pytorch operations and contain newly
implemented algorithms that are not available in standard Pytorch.

In this paper, we make a distinction between so-called native Pytorch
operations and Pytorch Geometric operations. Crucially, in cases where
analogous operations exist between the two (e.g. torch_scatter.scatter_add and
torch.Tensor.scatter_(op="add")), we provide side-by-side performance
benchmarks across several parameters such as data sparsity. We believe that
the resulting analysis will be useful to understand the optimizations provided
by the Pytorch Geometric library.

It is important to note that, while outside the scope of this paper,
alternatives to Pytorch Geometric exist, even within the Pytorch ecosystem.
For instance, Deep Graph Library (DGL) [5] provides native Python
implementations of message passing kernels, which, perhaps at the performance
costs of Pytorch/Pytorch Geometric’s native C++ implementations, provide
software platform agnostic support for message passing and other operations
important to GNNs. In the future, the performance of other GNN libraries such
as DGL can be compared to the existing Pytorch Geometric benchmarks found
below. Additionally, we plan to expand scope to other popular deep learning
frameworks such as TensorFlow.

###  I-C Hardware Support for Sparse Operations

A common problem facing many state-of-the-art deep learning systems (including
GNNs), is how to efficiently train and deploy models with millions, or even
billions, of learnable parameters [6]. As the size of neural networks continue
to increase, so does the importance of efficient and scalable hardware systems
to support them.

Over the past several years, NVIDIA GPUs have become a de-facto standard
hardware platform for training and deploying deep learning models at scale,
due primarily to their hardware parallelism and support for common deep
learning software tools, among other factors [7] [8] [9].

However, as models and their associated data continue to scale, hardware
bottlenecks have become increasingly problematic. In response, emerging
specialized hardware compute platforms have attempted solve these problems of
scalability and absolute performance. For example, the SambaNova RDU platform,
Graphcore IPU architecture, and Cerebras CS-2 platform, have all attempted
solve these efficiency bottlenecks by developing new dataflow architectures,
expanding on device memory and compute power, and/or providing explicit
hardware optimizations for deep learning [10] [11] [12] [13] [14]. Therefore,
rigorous operation-level performance benchmarks on NVIDIA systems can serve as
baseline performance references in order to facilitate future development of
such hardware systems.

###  I-D Profiling and Benchmark tools

Previous work has used several profiling and benchmarking frameworks in order
to evaluate both operations and applications (model) level performance of deep
learning models. For example, in the domain of GNNs, [15] developed a generic
benchmarking framework to facilitate applications level benchmarking of novel
model architectures. Other work used more ad hoc methods, or used native
Python tools such as the time.timeit module.

In this paper, we elect to use the native Pytorch profiler and native Pytorch
benchmark framework. This is due to several factors. Namely, since the scope
of this paper is limited to Pytorch and Pytorch Geometric, we are not limited
to generic profiling and benchmark tools. Additionally, the Pytorch profiler
supports automatic operations level profiles of end-to-end networks, which is
relevant to our native operation selection process (see III-A for more
details). The Pytorch benchmark framework extends the native Python timing
module to facilitate features such as automatic warmup runs, thread
synchronisation, and other details critical to equitable and accurate
benchmarking.

However, future work focusing on software frameworks other than Pytorch will
require a generic benchmark framework. Please refer to VI for more details.

##  II Related Work

In this section, we discuss prior work in profiling GNNs and benchmarking low-
level operations relevant to GNNs and other deep neural networks.

###  II-A Benchmarking Sparse Operations

Previous work has focused on benchmarking sparse operations at the kernel
level, benchmarking of specific novel sparse operations, and hardware
optimizations for sparse operators. For example, [16] demonstrates the
relative performance increases of NVIDIA’s cuSPARSE backend compared to the
standard cuBLAS, as well as proposes further kernel optimizations for certain
operations. Additionally, [17] proposed new software optimizations to improve
hardware efficiency of a particular operation, sparse-sparse matrix
multiplication. Other groups, such as [18], propose several hardware
architecture improvements to accelerate sparse operations, such as by unifying
memory management to support continguously represented sparse data. Mainstream
vendors such as NVIDIA have also proposed hardware improvements to their
systems in order facilitate more efficient sparse operations. However, to our
knowledge, no group has performed systematic benchmarks for low-level
operations critical to GNN performance at the systems level and within a
specific deep learning framework.

###  II-B GNN Benchmarking

Since the rise of popularity of GNNs, substantial research has focused on
creating effective and equitable benchmark paradigms. For instance, [15]
created one of the first generic and open source benchmark frameworks, wherein
researchers can easily add novel models and datasets and rigorously compare
results over end-to-end models, datasets, and hardware systems. Frameworks
such as these are critical for the development of more efficient GNN
architectures. However, such frameworks are generally concerned with end-to-
end model performance and stop short of providing operations level benchmarks
for GNNs.

Other previous work has focused on creating comprehensive and equitable
datasets in order to facilitate equitable end-to-end model performance.
Notably, [19] introduced Open Graph Benchmark, a set of diverse graph
benchmarks that can be used to assess model performance across several
metrics. These datasets are especially valuable in that they include both
realistic synthetic data that may be used for stress testing new models and
hardware systems as well as real world data that may be used for downstream
model performance analysis.

Perhaps most relevant to our work is the intial benchmarking conducted by [4]
, wherein Pytorch geometric implementations of several GNN architectures were
compared to existing implimentations, such as those of Deep Graph Library
(DGL) [5]. Additionally, [4] provided benchmark operations-level comparisons
of their scatter and gather operations. While they stopped short of
comprehensive operations level benchmarks, the results they did provide
enables us to robustly compare and resemblance check our results.

##  III Experimental Design

The benchmarking pipeline of this paper can be broken into two main sections.
First, common graph neural network (GNN) models are profiled in order to
determine key operations for GNNs (and sparse operations in general). We use
these results in order to collate these operations along with otheroperations
that are specifically designed to facilitate efficient sparse operations in
PyTorch Geometric.

Second, we systematically benchmark each of these operations on NVIDIA A100
GPUs. In the following subsections, we motivate and detail the procedure for
each of these steps.

###  III-A Model Profiling

The first motivation is to profile low-level operations that are integral to
modern GNN models. Thus, we select and robustly benchmark several end-to-end
GNN architectures. Relevant literature was reviewed in order to determine a
set of 8 models that we believe provide a prudent representation of the
current state-of-the-art in the field. Specifically, representations of
different scientific applications (molecules, social networks, traffic), graph
prediction task (node classification, node regression, graph regression, link
prediction, etc.), and architecture (graph attention, VAEs, spectral and
spatial methods) were considered. Additionally, models were screened by their
perceived impact on the field (e.g. citation count, use in industry, etc.) as
well as ease of implementation in the Pytorch Geometric framework. The final
models selected for profiling can be found in table I.

TABLE I: Selected GNN models for operations-level profiling

Model | Primary Domain
---|---
Graph Attention Network (GATv2) [20] | Proteins
Temporal Graph Network (TGN) [21] | Social Networks
Principal Neighbourhood Aggregation (PNA) [22] | Vision
Dynamic Neighborhood Aggregation (DNAConv) [23] | Several (representation learning)
Graph U-Net [24] | Several (representation learning)
SchNet [25] | Quantum Interactions
Directional Message Passing Network (DimeNet) [26] | Molecules

Each model was then profiled using the native Pytorch profiler, as discussed
in I-D. Each model was run on synthetic data of varying sparsity, first with
default model parameters. Then, the models were extended such that they
utilized full GPU memory. Specifically, models were extended by increasing
network parameters (iteratively deepening and widening each network). Where
extending networks lead to negligible performance increases, data batch sizes
were increased so as to reach full GPU utilization. Top operations were then
collated for each model and were used to select the native Pytorch operations
for benchmarking (IV-A2). Results from profiling these models are detailed in
IV.

###  III-B Operation Selection

In this section, we detail the process used to select relevant sparse
operations. As previously noted, Pytorch Geometric has four dependencies which
contain optimized operations for GNNs: torch_sparse, torch_scatter,
torch_cluster, and torch_spline_conv. In order to select operations for
benchmarking, we examined each of these libraries and selected fundamental
operations that were considered low-level and pertinent to GNNs. Crucially,
any available native Pytorch implementations of these operations were also
benchmarked in order to facilitate comparison between the two implementations.
Note that the torch_cluster library, which contains efficient implementations
of many high-level cluster operations, is not covered in this initial
benchmark experiment.

Additionally, the Pytorch operations which contributed the most to the overall
compute time across all models during label profiling (III-A) were evaluated
for criteria such as availability of sparse implementations for the CUDA
backend, relevance to other sparse operators, overlap with the operations
selected above, etc. Based on these constraints, five native Pytorch
operations were also selected for benchmarking.

The list of all operations selected for benchmarking can be found in IV-A1 and
IV-A2.

###  III-C Operation Benchmarking

In this section, we describe the experimental setup for operation
benchmarking. As discussed in II, experimental details of operations-level
benchmarking are rarely standardized, and to our knowledge, no standard
practice for this task exists. Thus, throughout this section, we justify
design decisions in-so-far as possible.

We test each operation at the full memory capacity of the machine in order to
facilitate fair comparison of operations across different hardware platforms.
We note that typical use cases of these operations may often use inputs far
smaller than those tested here. However, given the goal of many deep learning
accelerator hardware systems are to extend common GNN architectures,
benchmarking at GPU capacity facilitates comparison to hardware accelerator
systems that often have far higher memory capacity than standard GPUs. We
define this as $\geq 95\%$ memory capacity. Thus, for each pair of operation
and dimension (e.g. one dimensional tensors, two dimensional tensors, etc), we
select the input size so that the combination of input values with the highest
memory utilizes $\geq 95\%$ of the total memory capacity of the machine. All
other (less expensive) combinations of input values with the same dimension
for that operation will retain the same input size in order to allow for
robust comparisons of different operation parameters.

It is important to note that this later point implies that many of the input
combinations may yield peak memory loads of less than 95%. For example,
benchmarking a fully dense tensor will often require substantially more memory
overhead than a comparatively sparser one (even in cases where the they are
represented as strided tensors). However, in experiments where we benchmark an
operation over several sparsity values, we fix the input tensor to be of size
such that the most dense tensor utilizes full memory capacity. Thus, more
sparse tensors of the same size will not utilize full capacity. This enables
clear comparison of the effect of sparsity on the operation runtime.

After fixing the input size for different tensor shapes and varying tensor
sparsity, we also, where relevant, vary the other parameters of the operation.
However, in some cases, not all parameters are tested, especially in cases
where either relevant literature or our operation profiling indicated typical
use cases for an operation. We note that while more comprehensive data is
always beneficial, it is often not feasible to test and analyze all parameters
of every operation due to the polynomial increase in benchmark points required
for each new parameter.

Thus, the high-level process to determine the inputs to all operations is as
follows:

  1. 1.

Identify parameters to vary for the specific operation, always including
sparsity of the input tensor. Other parameters depend on the specific
operation (e.g. dimensionality of tensor, reduce operation, etc.).

  2. 2.

Given these combinations of input parameters, identify the combination (for
all input shapes) with highest peak memory capacity.

  3. 3.

Increase size of input tensors until this input combination reaches full
memory capacity ($\geq 95\%$).

  4. 4.

Set input tensor size for all other input combinations to the value above.

There are a number of additional assumptions and details associated with our
benchmarking method, such as data types, hardware details, etc. that are
summarized below:

  1. 1.

Pytorch version: 1.9.0

  2. 2.

CUDA version: 11.3

  3. 3.

Pytorch Geometric version: 2.0.4

  4. 4.

NVIDIA GPU details: NVIDIA A100-SXM4-40GB (40.536 DDR).

  5. 5.

All operations where benchmarked at least 14 times and all reported value are
the median value benchmarked.

  6. 6.

All operations benchmarked have an interquantile range of $\leq 0.00017$
seconds across single benchmark iterations.

  7. 7.

All tensors are assumed to be of type torch.float32 (data) or torch.int64
(indices) unless otherwise noted.

Additionally, we note that for this paper, we only benchmark the forward pass
of each operation, and do not include values for the backwards pass.

Finally, we define the term reduce factor (RF) to refer to the combination of
inputs needed in order for the ratio of the total elements of the input tensor
to the output tensor of an operation to be approximately the reduce factor.
For example, a reduce factor of 8 indicates that the output tensor of an
operation is about 8 times smaller than its input.

##  IV Results: Profiling and Operation Selection

###  IV-A Operation Selection

Below, we provide lists of the selected operations based on the method
described in the previous section.

####  IV-A1 Sparse Operations

The final operations selected for benchmarking from the torch_sparse,
torch_scatter, and torch_spline_conv libraries can be found in Table II below.

TABLE II: Selected Sparse Operations for Benchmarking

Operation | Library | Native Pytorch Counterpart | Notes
---|---|---|---
scatter_max | torch_scatter |  |
scatter_min | torch_scatter |  |
scatter_mean | torch_scatter |  |
scatter_add | torch_scatter | torch.Tensor.scatter_ (reduce=”add”) | Native impl. Inplace
|  | torch.Tensor.scatter_ (reduce=”multiply”) | Native impl. Inplace
coalesce | torch_sparse | torch.Tensor.coalesce |
transpose | torch_sparse | torch.Tensor.transpose | with COO formatted input
spmm | torch_sparse | torch.smm |
spspmm | torch_sparse | torch.sparse.mm |
spline_conv | torch_spline_conv |  |

####  IV-A2 Selected Pytorch/Aten Operations

The final operations selected for benchmarking from Pytorch can be found in
Table III below.

TABLE III: Selected Native Pytorch Operations for Benchmarking

Torch Level Operation Name | Analagous aten:: Operation | Notes
---|---|---
torch.addmm | aten::addmm |
torch.index_select | aten::index_select |
torch.sort | aten::sort |
torch.Tensor.index_add_ | aten::index_add_ | In-place tensor operation
torch.gather | aten::gather |

##  V Results: Operation Benchmarking

In the following subsections, the results of benchmarking the operations
listed in IV-A1 and IV-A2 are shown and analyzed. For each operation, or group
of operations, trends in data that were found to be particularly insightful
are shown and discussed. However, full benchmark numbers are omitted for
brevity and can instead be found at our project repository site
111https://github.com/ryienh/gnn-ops-benchmark. Additionally, visualizations
for some operations are omitted but can also be viewed in the project
repository site.

###  V-A Selected Native Pytorch Operations

This section explores the Pytorch operations selected as part of model
profiling in IV. These operations were found during profiling to be commonly
used in end-to-end GNN architectures but are not part of Pytorch Geometric or
its associated dependencies.

####  V-A1 torch.sort

Sorting operations have long been a focus of algorithms and
systems/performance research. The Pytorch sort operation sorts a single
strided input tensor. The documentation and source code for this operation can
be found here. Figure 1 provides a summary of the operation profile. Our
benchmarking experiments yield the following principal outcomes:

Figure 1:  Benchmarking results for the native torch.sort operation. This
operation is benchmarked for input dimensionalities between 1 and 3 and for
input sparsity values of 0, 0.5, 0.9, and 0.99. Additionally, the
implementations for stable and unstable sort are explicitly compared.

  1. (i)

In general, comparatively more sparse input tensors result in faster execution
times. This is probably due to the fact that sparser tensors have larger
chunks of zero elements that are already sorted.

  2. (ii)

The larger the dimensionality of the input tensor, the less the above effect
is observed, with the inverse correlation between input tensor sparsity and
operation runtime diminishing substantially. We hypothesize that this is due
to more intensive swap and copy operations, which bottleneck the algorithm and
thus reduce the significance of the input sparsity on the total runtime.

  3. (iii)

As expected, the stable sorting algorithm results in slower execution times
compared to its unstable counterpart.

  4. (iv)

As only densely formatted tensors are supported for this operation, the
sparsity of the input does not affect system memory.

  5. (v)

Sorting across higher dimensions of the input result in comparatively faster
execution times. This is consistent with the C/C++ style memory hierarchy.
However, in order to fully reconcile this phenomenon with the CUDA memory
hierarchy, further analysis focusing on information such as cache misses,
memory organization, etc. is necessary. Please see VI-C for more details.

Given the observations above, we believe that the primary bottlenecks that
should be targeted for future operation optimizations are related to memory
overheads. The implementation of a sort algorithm that operates on sparsely
formatted tensors will substantially improve memory and runtime overheads,
especially in multidimensional cases.

####  V-A2 torch.addmm

This operation combines a matrix multiplication with an add operation.
Specifically, the operation matrix multiplies two strided input tensors and
adds the result to a third strided input tensor. Thus, our benchmarking tests
different input combinations of 1 and 2 dimensions (using broadcasting where
applicable). The documentation and source code for this operation can be found
here. Figure 2 provides a summary of the operation profile. Our benchmarking
experiments yield the following principal outcomes:

  1. (i)

As all tensors are strided tensors, sparsity has no discernible effect on
runtime or memory performance.

  2. (ii)

As expected based on theoretical runtime analysis, performing matrix
multiplication on higher dimensional tensors result in substantially larger
runtimes. Varying the add tensor results in a far smaller, yet non-negligible
runtime increase.

The current multiply-and-add operation does not optimize for sparse tensor
input. Thus, implementation of a version of this operation for sparsely
formatted tensors will likely substantially improve memory overhead and
runtime performance for sparse tensor inputs.

Figure 2:  Benchmarking results for the native torch.addmm operation. This
operation is benchmarked for different combinations of input dimensionalities
between 1 and 2 and for input sparsity values of 0, 0.5, 0.9, and 0.99.

####  V-A3 torch.gather

At a high level, this operation gathers values of an input tensor by indexing
according to a given index tensor input. Our benchmarking test includes
different input combinations of 1, 2 and 3 dimensions and, in the latter two
cases, varies the dimension which is indexed. The documentation and source
code for this operation can be found here. Our benchmarking experiments yield
the following principal outcomes:

  1. (i)

Given the same total memory capacity, operations across a single dimension
result in comparatively larger runtimes vs. multidimensional tensors. This is
due to the fact that the absolute number of gather operations is larger when
gathering across a single dimension.

  2. (ii)

The above result extends to the total dimensionality of the input tensor more
generally. Specifically, when controlling for all other variables, the number
of dimensions of the input tensor is inversely proportional to the runtime.
This is a direct consequence of (i).

  3. (iii)

When gathering across a higher dimension, runtimes are comparatively longer.
We hypothesize that this is due to GPU memory storage hierarchy and access
factors. However, additional analysis in future work is required in order to
explore this further (See VI-C).

  4. (iv)

The index dimensionality parameter dominates the input dimensionality
parameter. Thus, low index dimensions with low dimensional tensor inputs
outperform high index dims with high dimensional tensor inputs.

  5. (v)

Due to strided tensor input formats, input sparsity has negligible effect on
total runtime and memory overhead.

Similar to other native Pytorch operations, we again observe negligible
effects of input tensor sparsity and runtime. Additionally, this operation
currently only supports strided tensor inputs, thus fixing the storage
capacity of the tensor regardless of input sparsity. Consequently, we believe
future implementation of a gather operation for sparsely formatted tensors
will greatly improve performance of this operation when used with sparse
inputs.

####  V-A4 torch.index_select

This operation uses an inputted index vector in order to index an inputted
value tensor along a given dimension. Our benchmarking test includes different
input combinations of 1, 2 and 3 dimensions and, in the latter two cases,
varies the dimension which is indexed. Additionally, we test several reduce
factor (RF) values, as defined in III-C. The documentation and source code for
this operation can be found here. Our benchmarking experiments yield the
following principal outcomes:

  1. (i)

Higher reduce factors lead to faster runtimes. This is due to the fact that
the number copy operations decreases linearly with the size of the output.
Additionally, the total operation memory overhead also decreases due to the
decreased output size of the operations at high reduce factor values.

  2. (ii)

Indexing at higher dimensions results in higher runtimes. Similar to V-A3,
this is likely due to memory hierarchy details of the GPU system, which is
outside the scope of this paper but should be explored further in future work.

  3. (iii)

Unlike, V-A3, reduce factor dominates index dimension. That is, high index
dimensions with high reduce factors outperform low index dimensions with low
reduce factors. This can be explained by comparing the indexing of the
torch.gather and torch.index_select operations: The former indexes across an
arbitrary number of index dimensions while the later uses the (1-dimensional)
index tensor to index the value tensor. Thus, the memory access order is
reversed.

  4. (iv)

Due to strided tensor input formats, sparsity has negligible effect on total
runtime and memory overhead.

We again observe negligible effects of input tensor sparsity and runtime due
to lack of support for sparsely formatted tensors and no explicit software
optimization for sparse tensors with strided format. Again, we believe the
implementation of a index_select operation for sparsely formatted tensors will
greatly improve performance of this operation when used with sparse inputs.

####  V-A5 torch.Tensor.index_add_

At a high level, this operation adds the value of an inputted source tensor to
the self tensor according to the indices given in the index input. Our
benchmarking test includes different input combinations of 1, 2 and 3
dimensions and, in the latter two cases, varies the dimension which is
indexed. The documentation and source code for this operation can be found
here. The associated visualizations for this operation are excluded from this
document for brevity but can be referenced in our project repository. Our
benchmarking experiments yield the following principal outcomes:

  1. (i)

Indexing at higher dimensions results in higher runtimes. This is consistent
with the memory hierarchy discussion outlined in V-A3 and V-A4.

  2. (ii)

There is no discernible pattern between dimensionality of input and the
runtime. This is a counterintuitive result that warrants future investigation
in future work. However, one hypothesis may be that discrepancies for total
element size of the input tensors may cause this issue. Specifically, the
memory overhead required for the multi-dimensional cases is such that there is
substantial $\geq 10\%$ discrepancies in total element size. Thus, future work
should normalize for total element size as well as investigate the memory sub-
system profiles further here.

  3. (iii)

Due to strided tensor input formats, sparsity has negligible effect on total
runtime and memory overhead.

As seen in the previous native operations, we again observe negligible effects
of input tensor sparsity and runtime due to lack of support for sparsely
formatted tensors and no explicit software optimization for sparse tensors
with strided format. Again, we believe the implementation of a index_select
operation for sparsely formatted tensors will greatly improve performance of
this operation when used with sparse inputs. Additionally, we encounter the
counterintuitve result wherein dimensionality of input has no correlation with
runtime performance. Future research is necessary in order to analyze and
verify this claim.

High-level Conclusions: The profiled native Pytorch operations only support
strided tensor inputs. With the exception of the sort operation (V-A1),
sparsity of the input tensors had negligible effects on runtime performance.
Given that these operations were shown in III-A to be a prominent part of
state-of-the-art GNNs, implementation of sparse equivalents of these
operations have the potential to improve end-to-end GNN performance,
especially when the input to these operations are sparse.

###  V-B Scatter operations

This section explores the benchmark results of scatter operations in the
torch_scatter package, one of the four dependencies of the Pytorch geometric
package. Scatter operations are a type of reduce operation that map an input
tensor to an output tensor by using an index tensor input. Where duplicate
indices exist, a reduce operation is used to select the final output value at
that index.

All scatter operations are benchmarked by varying the following parameters:
input dimension (1 and 2 dimensions), input sparsity, and reduce factor (as
defined in III-C). Where available, we provide direct comparisons to the
analogous native Pytorch operation, torch.Tensor.scatter_. We also note that
unlike the torch_scatter operations, torch.Tensor.scatter_ operations are in-
place.

The documentation and source code for these operations can be found at the
following hyperlinks: torch_scatter and torch.Tensor.scatter_.

####  V-B1 scatter_min and scatter_max

We first evaluate the results of benchmarking the scatter_min and scatter_max
operations. The scatter_min is visualized in Figure 3.

  1. (i)

In the one-dimensional case, higher reduce factors lead to comparatively
faster runtimes and smaller total operation memory overhead. Similar to V-A4,
this is due to the fact that the number copy operations decreases linearly
with the size of the output.

  2. (ii)

In the two dimensional case, higher reduce factors lead to comparatively
slower runtimes. This is may be due to the fact that the reduce operations are
calculated across two dimensions and are thus non-contiguous. However,
additional examination of the GPU memory system is necessary in order to
confirm this hypothesis.

  3. (iii)

Runtime and sparsity have a direct relationship, which is unexpected. Change
in runtime with respect to sparsity is most pronounced in low dimensional
cases. This may be due to the fact that higher dimensional tensors require far
more copy operations per step (one must copy entire tensor at $dim\neq
idx\\_dim$ to output). Thus, this bottlenecks the algorithm.

  4. (iv)

As expected, scatter_min and scatter_max have similar trends in runtime and
memory peaks due to the reduce operations being nearly identical in
computation time.

The observations above imply that sparse inputs lead to comparatively inferior
runtime speeds, perhaps due to memory non-contiguity and irregular accesses.
Additionally, the effect of the reduce factor on runtime varies depending on
the dimensionality of the input. Thus, we believe that memory hierarchy
optimizations may be the most promising course of action to improving
performance the performance of these operations in end-to-end GNN models.

Figure 3:  Benchmarking results for the scatter_min operation. This operation
is benchmarked for inputs of both 1 and 2 dimensions, reduce factors of 1, 2,
4, and 8, and input sparsity values of 0, 0.5, 0.9, and 0.99.

####  V-B2 scatter_mean and scatter_add

We now evaluate the results of benchmarking the scatter_mean and scatter_add
operations. Additionally, scatter_mean is visualized in Figure 4. We note that
the analogous native Pytorch implementation to scatter_add is
torch.Tensor.scatter_ (reduce="add"), and both operations’ results are
included in our analysis.

  1. (i)

Both scatter_mean and scatter_add have negligible changes in runtime as input
sparsity is varied. We hypothesize that this is due to mean and add operations
requiring more computation compared to min and max. Additionally, analysis of
the source code indicates that optimization for 0 valued elements is present
in the latter but not the former.

  2. (ii)

We also compare scatter_add with the native Pytorch operation
torch.Tensor.scatter_ (reduce="add"). Performance of the native operation is
comparatively faster than than scatter_add. This is likely largely due to the
in-place nature of the native operation.

Unlike with the min and max scatter operations, varying the sparsity of the
input has only negligible effects on overall runtime. Thus, we believe the
biggest computational bottlenecks for these operations are either the reduce
operations themselves or the memory hierarchy of the GPU.

Figure 4:  Benchmarking results for the scatter_mean operation. This operation
is benchmarked for inputs of both 1 and 2 dimensions, reduce factors of 1, 2,
4, and 8, and input sparsity values of 0, 0.5, 0.9, and 0.99.

####  V-B3 torch.Tensor.scatter_ (reduce="multiply")

We also benchmark the native torch.Tensor.scatter_ (reduce="multiply"), which
has no torch_scatter analogue. We note that as with the earlier add and mean
scatter operations, varying the sparsity of torch.Tensor.scatter_
(reduce="multiply") has a negligible effect on the total runtime and that the
effect of the reduce factor depends on the dimensionality of the input.

High-level Conclusions The chosen reduction operation (e.g. min, max, add,
etc.) of benchmarked scatter operations significantly influence the effect of
other factors, such as sparsity and input dimensionality of the input, on the
overall runtime. Examination of source code indicates significantly different
implementation decisions based on the reduction operation selected. Thus,
simplification to a single efficient implementation may improve performance of
many of these operations.

Additionally, as described in the above sections, significant variance in
runtime is observed across different input dimensionalities (despite
approximately constant total number of elements and reduce operations
required). Thus, these operations may benefit from memory hierarchy designs
that are developed for sparse inputs.

###  V-C Sparse operations

This section explores the benchmark results of miscellaneous sparse operations
in the torch_sparse package, one of the dependencies of the Pytorch geometric
package. Specifically, we examine sparse-dense and sparse-sparse matrix
multiplication, sparse matrix transpose, and sparse coalesce operations. All
of the operations in this section contain analogous native Pytorch operations,
which are included in the analysis. We note that all native Pytorch operations
use sparse (COO) formatted tensors.

####  V-C1 sparse_spmm and sparse_spspmm

We first evaluate the results of benchmarking the sparse_spmm and
sparse_spspmm operations, which correspond to sparse-dense and sparse-sparse
matrix multiplication respectively. Both operations are analogous to
torch.sparse.mm (as this operation can operate on both COO and strided tensor
inputs). Thus, we compare sparse_spmm to torch.sparse.mm with a single COO and
single strided tensor input and sparse_spspmm to torch.sparse.mm with two COO
formatted inputs. These results are visualized in Figures 5 and 6,
respectively. The documentation and source code for these operations can be
found here, and for the native operation here.

Figure 5:  Benchmarking results for the sparse_spmm operation. This operation
is benchmarked for combinations of inputs of 1 and 2 dimensions and input
sparsity values of 0, 0.5, 0.9, and 0.99. Solid lines and small markers denote
PyTorch Geometric operations, whereas dashed lines and large markers denote
PyTorch operations.

  1. (i)

sparse_spmm underperforms its native equivalent in all but the 1 dimensional
input case.

  2. (ii)

In general, sparse_spspmm also underperforms its native equivalent, across all
sparsity levels.

  3. (iii)

The sparsity of the first tensor (for sparse-dense multiplication) and the
overall sparsity (for sparse-sparse multiplication) is inversely proportional
to both runtime and memory requirements.

  4. (iv)

This effect is most pronounced at high sparsity levels (e.g. between 90% and
99% sparsity). Thus, future work can investigate performance differences
across sparsity levels close to 100%.

Figure 6:  Benchmarking results for the sparse_spspmm operation. This
operation is benchmarked for combinations of inputs of 1 and 2 dimensions and
input sparsity values of 0, 0.5, 0.9, and 0.99. Small markers denote PyTorch
Geometric operations, large markers denote PyTorch operations.

The above results imply that the native Pytorch torch.sparse.mm operation
often outperforms its Pytorch Geometric counterpart. Thus, the runtime
performance of end-to-end GNN models on GPUs will likely improve with either
further optimization of the Pytorch Geometric operations or use of the native
Pytorch operations.

####  V-C2 sparse_transpose

We now evaluate the sparse matrix transpose operation sparse_transpose and its
native Pytorch analogue torch.transpose. The documentation and source code for
this operation and its analogous Pytorch native equivalent can be found here
and here, respectively. We highlight the following observations:

  1. (i)

The Pytorch Geometric transpose operation under-performs its native
counterpart. This effect is more pronounced at comparatively lower input
sparsity values.

  2. (ii)

For both operations, runtime increases logarithmically with the input
sparsity.

  3. (iii)

The Pytorch Geometric operation runtime performance suffers significantly in
the two-dimensional input case. For example, for a two-dimensional input size
and sparsity of 50%, the Pytorch Geometric operation is on the order of 200
times slower than its native Pytorch counterpart. We hypothesize that this may
be due to lack of memory optimization for copy operations for the Pytorch
Geometric operation.

Based on these results, we conclude that, as with the sparse matrix
multiplication operations, that the use of native Pytorch operations will
likely improve end-to-end GNN performance.

####  V-C3 sparse_coalesce

Next, we evaluate the coalesce operation sparse_coalesce and its native
Pytorch analogue torch.coalesce. These results are visualized in Figure 7. The
documentation and source code for this operation and its analogous Pytorch
native equivalent can be found here and here, respectively. We highlight the
following observations:

Figure 7:  Benchmarking results for the sparse_coalesce operation. This
operation is benchmarked for inputs of 1 and 2 dimensions, reduce factors of
1, 2, 4, and 8, and input sparsity values of 0, 0.5, 0.9, and 0.99. Solid
lines and small markers denote PyTorch Geometric operations, whereas dashed
lines and large markers denote PyTorch operations.

  1. (i)

For both operations, runtime decreases with the input sparsity. This effect is
especially pronounced between 90% and 99% input sparsity levels.

  2. (ii)

For all input shapes, a higher reduce factor leads to comparatively higher
runtimes. This is due to the fact that a higher reduce factor requires a
comparatively higher number of add operations in order to coalesce the tensor.

  3. (iii)

For all but the trivial case where the reduce factor is 1, the native Pytorch
Geometric operation under-performs its native Pytorch equivalent. This result
holds across all input tensor sparsity levels.

We again conclude that the native Pytorch implementation is competitive with
the Pytorch Geometric implementation on GPU and should be considered in order
to improve GNN runtime performance.

High-level Conclusions: Based on the torch_sparse operations benchmarking, we
note that, unlike in previous operation groups, a key bottleneck is input
tensor sparsity. Thus, further software and hardware optimization should focus
on efficiency of especially sparse (and specifically COO formatted) tensors.
Additionally, we find that for all operations benchmarked, native Pytorch
operations are either competitive with, or outperform, the Pytorch Geometric
operations. Thus, the use of these operations should be considered as part of
a larger GNN implementation and future work should focus on optimization of
the Pytorch Geometric operations.

###  V-D Spline Convolution Operation

This section focuses on benchmarking the Pytorch Geometric implementation of
the spline convolution operator. The spline convolution operator is a
generalization of deep neural networks that operates on graph structured data,
and, by utilizing B-splines, is able to decouple kernel size from computation
time [27]. The associated Pytorch Geometric operation is
torch_spline_conv.spline_conv and its documentation can be found here.

We benchmark this operator over input sparsity and kernel size parameters. We
note that, unlike the preceding operations, sparsity is presented in the form
of average node degree which is inversely proportional to sparsity.
Additionally, due to differences in total node size in our stochastically
generated synthetic dataset, we present our results as GPU time normalized by
the total node size of the input. Our key observations can be found below:

Figure 8:  Benchmarking results for the spline_conv operation. This operation
is benchmarked for kernel sizes of 1, 3, 5 and 7, and input sparsity values of
0, 0.5, 0.9, and 0.99.

  1. (i)

Total runtime is directly proportional to the average node degree and thus
inversely proportional to the input sparsity.

  2. (ii)

Larger convolution kernels result in comparatively longer absolute runtimes,
especially at higher levels of input data sparsity.

While the above observations are fairly intuitive, we hope that the absolute
benchmark numbers provided will help facilitate further hardware and software
based optimization of this operation.

##  VI Future Work

While we believe that the operation profiling and benchmark framework proposed
in this paper, along with the results and associated analysis and insights, is
a valuable contribution, it is not comprehensive over all operations relevant
to GNNs. Rather, we believe it lays the groundwork for future study of
operations level performance across different software frameworks, hardware
platforms, model architectures, etc. In this section, we outline areas of
improvement for future work in the area.

###  VI-A Operation Benchmarks in Training Regime

This paper focuses solely on the forward pass of each operation during
benchmarking. Thus, performance is limited to the inference stage of GNNs.
Future work should build on our framework in order to facilitate benchmarking
of the backwards pass implementations of each operation as well. We believe
this should be a relatively simple extension of the current implementation of
the benchmarking framework.

###  VI-B Expanding Operation Coverage

Future work can expand on the results and analysis found here for other sparse
operations. For example, the torch_cluster library and other top operations
profiled in III-A may serve as starting points for further research in this
direction.

###  VI-C Additional Benchmark Metrics

Throughout this paper, several hypotheses were proposed to explain empirical
results of each operation. Thorough investigation of these hypotheses will
require benchmarking with additional systems metrics, such as cache misses,
data representation in memory, and specific hardware bottlenecks in each
operation.

###  VI-D Extending Comparisons to Other Hardware Platforms

We hope that the GPU results described in this paper serve as baseline
performances for those developing these operations on specialized hardware.
Additionally, our analysis may also prove valuable to those developing
optimized hardware architectures to support sparse operations. Future work may
focus on benchmarking these operations on existing hardware platforms in order
to evaluate their performances compared to a single NVIDIA A100 machine.

## Conclusion

In conclusion, we provide the following high-level analysis of our operation
benchmarks, summarizing the operation level observations of section (V):

  1. 1.

In several cases, such as in V-A and V-B, memory hierarchy order dominates
sparsity. Thus, in these cases, bottlenecks identified for optimization should
perhaps focus on memory, rather than solely input data sparsity.

  2. 2.

Native Pytorch operations should be considered competitive, often matching or
outperforming their Pytorch Geometric counterparts. This is especially true
when data is not very sparse, and holds even when the native Pytorch
equivalent is not an in-place operation. This performance gap is especially
pronounced for sparse operations such as transpose and sparse multiply, where
native operations support COO formatted sparse tensor representations.

  3. 3.

While some ops with strided inputs are optimized for comparatively sparser
inputs (e.g. torch.sort), many (such torch.addmm) have no such optimizations.
This can cause both memory and runtime bottlenecks which can significantly
affect end-to-end GNN model performance. Thus, GNNs will likely benefit from
sparse implementations of the benchmarked native operations.

We plan to expand our efforts to cover additional accelerator architectures
and model architectures in the future.

## Acknowledgements

This research used resources of the Argonne Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC02-06CH11357.

## References

  * [1] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu.  A comprehensive survey on graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 32(1):4–24, 2021.
  * [2] Michael M. Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst.  Geometric deep learning: Going beyond euclidean data.  IEEE Signal Processing Magazine, 34(4):18–42, 2017.
  * [3] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32, 2019.
  * [4] Matthias Fey and Jan Eric Lenssen.  Fast graph representation learning with pytorch geometric.  arXiv preprint arXiv:1903.02428, 2019.
  * [5] Minjie Yu Wang.  Deep graph library: Towards efficient and scalable deep learning on graphs.  In ICLR workshop on representation learning on graphs and manifolds, 2019.
  * [6] Maurizio Capra, Beatrice Bussolino, Alberto Marchisio, Guido Masera, Maurizio Martina, and Muhammad Shafique.  Hardware and software optimizations for accelerating deep neural networks: Survey of current trends, challenges, and the road ahead.  IEEE Access, 8:225134–225180, 2020.
  * [7] Aniruddha Parvat, Jai Chavan, Siddhesh Kadam, Souradeep Dev, and Vidhi Pathak.  A survey of deep-learning frameworks.  In 2017 International Conference on Inventive Systems and Control (ICISC), pages 1–7. IEEE, 2017.
  * [8] Giang Nguyen, Stefan Dlugolinsky, Martin Bobák, Viet Tran, Alvaro Lopez Garcia, Ignacio Heredia, Peter Malík, and Ladislav Hluchỳ.  Machine learning and deep learning frameworks and libraries for large-scale data mining: a survey.  Artificial Intelligence Review, 52(1):77–124, 2019.
  * [9] Sparsh Mittal and Shraiysh Vaishay.  A survey of techniques for optimizing deep learning on gpus.  Journal of Systems Architecture, 99:101635, 2019.
  * [10] Raghu Prabhakar and Sumti Jairath.  Sambanova sn10 rdu: Accelerating software 2.0 with dataflow.  In 2021 IEEE Hot Chips 33 Symposium (HCS), pages 1–37. IEEE, 2021\.
  * [11] M. Emani, V. Vishwanath, C. Adams, M. E. Papka, R. Stevens, L. Florescu, S. Jairath, W. Liu, T. Nama, and A. Sujeeth.  Accelerating scientific applications with sambanova reconfigurable dataflow architecture.  Computing in Science & Engineering, 23(02):114–119, mar 2021.
  * [12] Thorben Louw and Simon McIntosh-Smith.  Using the graphcore ipu for traditional hpc applications.  In 3rd Workshop on Accelerated Machine Learning (AccML), 2021.
  * [13] Lakshan Ram Madhan Mohan, Alexander Marshall, Samuel Maddrell-Mander, Daniel O’Hanlon, Konstantinos Petridis, Jonas Rademacker, Victoria Rege, and Alexander Titterton.  Studying the potential of graphcore ipus for applications in particle physics.  arXiv preprint arXiv:2008.09210, 2020.
  * [14] Asif Mirza, Shadi Manafi Avari, Ebadollah Taheri, Sudeep Pasricha, and Mahdi Nikdast.  Opportunities for cross-layer design in high-performance computing systems with integrated silicon photonic networks.  In 2020 Design, Automation Test in Europe Conference Exhibition (DATE), pages 1622–1627, 2020.
  * [15] Vijay Prakash Dwivedi, Chaitanya K Joshi, Thomas Laurent, Yoshua Bengio, and Xavier Bresson.  Benchmarking graph neural networks.  arXiv preprint arXiv:2003.00982, 2020.
  * [16] Trevor Gale, Matei Zaharia, Cliff Young, and Erich Elsen.  Sparse gpu kernels for deep learning.  In SC20: International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–14. IEEE, 2020.
  * [17] Guyue Huang, Guohao Dai, Yu Wang, and Huazhong Yang.  Ge-spmm: General-purpose sparse matrix-matrix multiplication on gpus for graph neural networks.  In SC20: International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–12. IEEE, 2020.
  * [18] Shijin Zhang, Zidong Du, Lei Zhang, Huiying Lan, Shaoli Liu, Ling Li, Qi Guo, Tianshi Chen, and Yunji Chen.  Cambricon-x: An accelerator for sparse neural networks.  In 2016 49th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO), pages 1–12. IEEE, 2016.
  * [19] Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  Advances in neural information processing systems, 33:22118–22133, 2020.
  * [20] Shaked Brody, Uri Alon, and Eran Yahav.  How attentive are graph attention networks?  In International Conference on Learning Representations, 2021.
  * [21] Emanuele Rossi, Ben Chamberlain, Fabrizio Frasca, Davide Eynard, Federico Monti, and Michael Bronstein.  Temporal graph networks for deep learning on dynamic graphs.  arXiv preprint arXiv:2006.10637, 2020.
  * [22] Gabriele Corso, Luca Cavalleri, Dominique Beaini, Pietro Liò, and Petar Veličković.  Principal neighbourhood aggregation for graph nets.  Advances in Neural Information Processing Systems, 33:13260–13271, 2020.
  * [23] Matthias Fey.  Just jump: Dynamic neighborhood aggregation in graph neural networks, 2019\.
  * [24] Hongyang Gao and Shuiwang Ji.  Graph u-nets.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 2083–2092. PMLR, 09–15 Jun 2019.
  * [25] Kristof Schütt, Pieter-Jan Kindermans, Huziel Enoc Sauceda Felix, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller.  Schnet: A continuous-filter convolutional neural network for modeling quantum interactions.  Advances in neural information processing systems, 30, 2017.
  * [26] Johannes Klicpera, Janek Groß, and Stephan Günnemann.  Directional message passing for molecular graphs.  In International Conference on Learning Representations, 2020.
  * [27] Matthias Fey, Jan Eric Lenssen, Frank Weichert, and Heinrich Müller.  Splinecnn: Fast geometric deep learning with continuous b-spline kernels.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 869–877, 2018.