Accelerating Graph Neural Networks with a Novel Matrix Compression Format
João N. F. Alves Samir Moustafa Siegfried Benkner Alexandre P. Francisco Wilfried N. Gansterer Luís M. S. Russo

[
    [
    

=====



§ ABSTRACT 
  The inference and training stages of Graph Neural Networks (GNNs) are often dominated by the time required to compute a long sequence of matrix multiplications between the sparse graph adjacency matrix and its embedding. To accelerate these stages, we first propose the Compressed Binary Matrix (CBM) storage format to succinctly represent the binary adjacency matrix of an unweighted graph. Then, we show how to generalize this representation to normalized adjacency matrices of unweighted graphs which arise in the context of GNNs. Finally, we develop efficient matrix multiplication kernels based on this compressed representation. The matrix multiplication kernels proposed in this work never require more scalar operations than classic sparse matrix multiplication algorithms. Experimental evaluation shows that the matrix multiplication strategies proposed outperform the current state-of-the-art implementations provided by Intel MKL, achieving speedups close to 5×. Furthermore, our optimized matrix-multiplication strategies accelerated the inference time of a GNN by up to 3×.









§ INTRODUCTION


Graph Neural Networks (GNNs) are the preferred tool to learn from
graph-structured data and thus are considered key for future AI applications in
domains like social network analysis, natural language processing, biology,
physics, and many others <cit.>. The training and inference time of different GNN architectures is dominated by a long sequence of matrix products. This is particularly evident in GNNs that resort to Message Passing Layers (MPLs), where in each hidden layer the nodes of the graph aggregate the embedding of neighboring nodes and adjust their own embedding based on the information collected. 
In some variants of GNNs, such as the widely used Graph
Convolutional Networks (GCNs) <cit.>, the message produced in each layer is
essentially the product of the adjacency matrix of the underlying graph and its current embedding. 

For illustration purposes, consider a two-layer GCN. To propagate and combine node features across the graph, this network must compute the following operations once per inference and once per training epoch <cit.>:

    𝐀̂ σ(𝐀̂𝐗𝐖^0) 𝐖^1,
 where 𝐀̂ represents the normalized adjacency matrix of the 
graph, such that 𝐀̂ = 𝐃^-1/2 (𝐀 + 𝐈) 𝐃^-1/2, 𝐃 is
the degree diagonal matrix of the graph, σ denotes an element-wise activation function, 𝐗 is the matrix of node features, and 𝐖^0 and 𝐖^1 are learnable dense matrices for the first and second layers <cit.>. It is important to note that in real scenarios the adjacency matrix of the graph is typically much larger than the remaining operand matrices of Equation <ref>. Therefore, matrix products involving this matrix represent most of the computational burden of training and inferring GNNs. Graphs arising in the context of GNNs are often extremely sparse. Popular GNNs frameworks, such as PyTorch <cit.>, leverage this sparsity to accelerate both training and inference. This is achieved by
representing the adjacency matrix of the graph, or its normalized form, in standard sparse matrix formats that only consider the nonzero elements of a sparse matrix. Sparse matrix formats, like COO or CSR, enable faster sparse-dense matrix multiplication kernels (SpMMs) which are known to require a number of scalar operations that is proportional to the number of nonzero elements in the sparse matrix.

In this paper we present a new matrix compression format called Compressed Binary Matrix (CBM) to accelerate widely-used GNN architectures (e.g., GCN <cit.>, GraphSage<cit.>, and GIN <cit.>)  that spend most of their processing time on matrix products involving the adjacency matrix of the underlying graph or its normalized, when these GNNs learn on unweighted graphs. Our format exploits the fact that the adjacency matrix of an unweighted graph is binary, and therefore it can be compressed beyond what is possible with sparsity alone, while simultaneously reducing the number of scalar operation required to multiply our compressed representation of the adjacency matrix and another dense real-valued matrix.
The key advantage of the CBM format is that it only represents the differences (deltas) of a row with respect to another similar row of the same matrix, which for many use cases tends to be significantly smaller than the number of non-zeros represented by standard compression formats.



 §.§ Previous Works


The matrix-matrix and matrix-vector products have been extensively studied <cit.>.
A particular case is the product of a binary sparse matrix by a real-valued vector or (dense) matrix, where the efficient representation of the binary matrix can be exploited to improve both the memory footprint and the operation running time.
Although this was not expected in some preliminary studies <cit.>, and impossibility results exist for more complex compression schemes <cit.>, it works for some representational compression schemes.

One of such schemes is the Single Tree Adjacency Forest (STAF) proposed by Nishino et al. <cit.>. The STAF representation of a binary matrix is obtained by reversing and inserting the adjacency list of each row of a binary matrix into a trie data-structure, meaning that suffixes that are common to more than one adjacency list are represented exactly once. The authors have shown that the this data-structure can be built in linear time with respect to the number of nonzero elements in the input matrix.
STAFs enable fast matrix-matrix products by traversing the trie in root-to-leaf order, while accumulating partial sums that are common to different rows of the result matrix.
The number of operations required to multiply a binary matrix represented by a STAF and a real-valued vector is proportional to the size of the trie, which is also upper-bounded by the number of nonzero elements in the input matrix.

STAFs doe not exploit, however, row-wise similarities outside of common row suffixes. To address this issue the authors proposed splitting the input matrix into sets of columns and create a STAF per set. This optimized version achieves a significant speedup and memory footprint reduction compared to CSR and the sparse matrix multiply kernel offered by the Eigen library. 

Francisco et al. <cit.> explored how succinct representations for binary matrices and graphs could be exploited to speedup binary matrix multiplication, namely the Webgraph representation by Boldi and Vigna <cit.> and the Biclique Extraction (BE) based representation by Hernandez and Navarro <cit.>.
Both representations allow to reduce the memory footprint of binary matrices and accelerate the product of compressed binary matrices and real-valued vectors.
The Webgraph representation exploits the similarity among rows, as well as clustering effects present in real-world graphs and matrices, relying on gap compression, referentiation, intervalisation and ζ codes.
The BE representation exploits also rows similarities and clustering effects, extracting maximal bicliques and replacing them with differential compressed sets of nodes that share adjacencies.


The key observations are that, in both cases, we can reuse partial results from previous computations.
Taking an implementation of Pagerank using classic adjacency lists, authors achieved significant and similar speedups and memory footprint reductions, showing that the product computation can be reduced in time proportional to the compressed matrix size. 
The Webgraph and BE representations require however non-trivial preprocessing steps, Webgraph benefits from a suitable vertex reordering, obtained through graph clustering methods, and BE requires finding maximal bicliques, an NP-hard problem in general.
These steps might take longer to execute than computing the SpMM we are trying to optimize.
On the other hand the experimental evaluations considered naive implementations of standard representations, possibly leading to unfair comparisons in what concerns state of the art linear algebra libraries.

Elgohary et al. <cit.> also addressed the problem that large-scale machine learning algorithms are often iterative, using repeated read-only data access and I/O-bound matrix-vector multiplications, introducing Compressed Linear Algebra (CLA) for lossless matrix compression.
CLA also executes linear algebra operations directly on the compressed representations.
However, it is not focused on binary matrices and, although it achieves a performance close to the uncompressed case, it only presents performance gains when data does not fit into memory.

Our work is somewhat related to the work by Björklund and Lingas <cit.>.
They also consider a weighted graph on the rows of a binary matrix where the weight of an edge between two rows is equal to its Hamming distance, and then they rely on a minimum spanning tree of that graph to differentially compress the rows of a binary matrix with respect to one another, thus accelerating the product computation.
They consider however only a single product of two binary matrices, and do not consider the overhead imposed by operating on a compressed representation of a binary matrix, as their results are purely theoretical.



 §.§ Our Contributions

In this paper, we make the following contributions: (1)  we present the Compressed Binary Matrix (CBM) format, an efficient compression scheme for binary matrices, that can reduce the memory footprint of unweighted graphs that arise in the context of GNNs; (2) we introduce a new algorithm to significantly accelerate the sequence of matrix products between the (potentially normalized) adjacency matrix of the graph, represented in CBM format, and its embedding; (3) we prove that even in the worst-case scenario, where compression
is not possible, the number of scalar operations required to multiply a matrix represented in CBM format does not exceed the number of scalar operations required to multiply the same matrix using classic sparse storage formats; and (4) we have implemented the CBM format and corresponding matrix multiplication kernels such that they can be used together with state-of-the-art Deep Learning framework, such as PyTorch. Furthermore, experimental evaluation using real-world datasets demonstrates the effectiveness of our approach. Our method is nearly 5× faster than state-of-the-art SpMM implementations in sequential and parallel environments, subsequently shortening the inference time of a 2-layer GCN by more than 3×. Our implementation will be made available at <https://github.com/cbm4scale>.

As previously stated, the CBM format is akin to the work of Björklund and Lingas <cit.>. Nevertheless, it distinguishes itself by being specifically designed to accelerate the product between a single sparse binary matrix and a set of real-valued matrices, thus the binary matrix only needs to be compressed once[Ideally, the adjacency matrix of any unweighted graph would be provided in CBM, avoiding any compression overhead.]. Furthermore, our format resorts to a Minimum Cost Arborescence (MCA), which allows us to ignore compression opportunities that do not lead to improvements in memory-footprint and matrix multiplication. The CBM format also overcomes limitations present in STAF and BE, since our format exploits compression opportunities along complete rows of the matrix (by default), is constructed in polynomial time, and does not allocate additional memory to store partial results. 



§ COMPRESSED BINARY MATRIX FORMAT

Let 𝐀∈{0,1}^m × n be a binary matrix, where a⃗_i ∈{0,1}^n represents the i-th row of 𝐀 for i=1, …, m. Also assume that a⃗_i is represented with an adjacency list with the column indices of the nonzero elements of a⃗_i.
The Compressed Binary Matrix (CBM) format resorts to differential compression to represent the rows of a binary matrix 𝐀. This is, if 𝐀 is represented in CBM format, then any row a⃗_x can be represented by another row a⃗_y and two lists of deltas that indicate which elements must be included (Δ^+_x,y), or removed (Δ^-_x,y), from the adjacency list of a⃗_y to obtain a⃗_x:

    a⃗_x = (a⃗_y ∪Δ^+_x,y )∖Δ^-_x,y,  where Δ^+_x,y=a⃗_x ∖a⃗_y , and Δ^-_x,y=a⃗_y ∖a⃗_x.


Assuming a⃗_y is present in memory, Equation <ref> suggests that the memory required to represent a⃗_x is proportional to the number of deltas between a⃗_x and a⃗_y. If these rows are similar, then it is likely that the number of deltas is smaller than the number of nonzero elements of a⃗_x. In that case, it would be more memory efficient to represent a⃗_x with respect to a⃗_y than through its adjacency list.
Therefore, to further reduce the memory footprint of matrix A, the compression algorithm that builds the CBM format must find a suitable chain of compression to represent all rows of 𝐀. This is, for each row a⃗_x, identify another similar row a⃗_y that characterizes the former, such that: (1) the numbers of deltas required to represent each row a⃗_x is minimized subjected to a⃗_y, and (2) the number of deltas required to represent a_x is guaranteed to be less than, or equal to, the number of nonzero elements in a⃗_x.



  
Minimizing the number of deltas. To address point (1), the CBM format must first measure the number of deltas required to convert each row a⃗_y into all other rows a⃗_x of 𝐀, i.e., measure the Hamming distance for each pair of matrix rows.
This step provides a global view of the dissimilarity between the rows of the matrix 𝐀, and it can be modeled as a fully-connected and undirected distance graph G. This graph has m vertices, where each vertex represents a unique row of the matrix, and the weight of each edge (y,x) corresponds to the number of deltas required to represent a⃗_x with respect to a⃗_y.
To reduce the number of deltas required to compress the rows of 𝐀 we can find a Minimal Spanning Tree (MST) of G, which by definition spans G with the minimum sum of edge weights possible. Naturally, any MST of the distance graph, rooted in vertex x, defines a chain of compression with as many deltas as the weight of this tree plus the number of nonzero elements of a⃗_x. Thus, any MST rooted in the vertex corresponding to the row with the fewest nonzero elements, defines a chain of compression that satisfies point (1).



  
Worst-case guarantees. Note that the chain of compression obtained by finding an MST of G does not satisfy point (2), because the weight of the lightest incoming edge of any vertex x can be greater than the number of nonzero elements in a⃗_x. In such cases, representing a⃗_x with an adjacency list is clearly more memory efficient. To avoid this issue, we extended the distance graph G with a virtual vertex 0 which is connected to all other vertices of the graph. This virtual vertex represents a null row-vector a⃗_0 ∈{0}^n, which ensures that the weight of each edge (0,x) is equal to the number of nonzero elements in a⃗_x. The inclusion of virtual vertex 0 in the distance graph G ensures that the issue described above cannot occur, since the lightest incoming edge of any vertex x is now at most as heavy as the number of nonzero elements in a⃗_x. Therefore, any chain of compression characterized by an MST of G, rooted on vertex 0, is guaranteed to satisfy points (1) and (2). If we use this chain of compression to represent a matrix in CBM format, then following property will be observed:


    The number of deltas required to represent any matrix 𝐀∈{0,1}^m × n in Compressed Binary Matrix (CBM) format is never greater than the number of nonzero elements in matrix 𝐀.


To complete the construction, we simply traverse the compression chain above in topological order, and for every edge (y,x) visited, we compute the lists of positive and negative deltas required to convert row a⃗_y into a⃗_x.



 §.§ Time and Space Analysis



Any matrix 𝐀∈{0,1}^m × n can be represented in CBM format in O((m+1) 𝐧𝐧𝐳(𝐀) + m^2 log m) time. 

    The construction of the extended distance graph G for 𝐀 requires the computation of m(m+1)/2 Hamming distances between all possible row pairs (a⃗_x, a⃗_y). The Hamming distance of each pair of rows can be reduced to the intersection of their adjacency lists, computed in O(𝐧𝐧𝐳(a⃗_x) + 𝐧𝐧𝐳(a⃗_y)) time. Hence, the time to compute all m(m+1)/2 Hamming distances is upper-bounded by
    
    ∑_x=0^m ∑_y=0^m (𝐧𝐧𝐳(a⃗_x) + 𝐧𝐧𝐳(a⃗_y)) = (m+1) 𝐧𝐧𝐳(𝐀).

Additionally, well-known MST algorithms, such as Prim or Kruskal, are known to find an MST in O(E log V) time, where E and V denote the number of edges and vertices in the graph. Since the extended distance graph G contains m(m+1)/2 edges and m+1 vertices, finding an MST of G requires O((m+1)^2 log (m+1)) time, and therefore, representing matrix 𝐀 takes time

    O((m+1) 𝐧𝐧𝐳(𝐀) + (m+1)^2 log (m+1))).


Finally, a single list of deltas can also be obtained from the intersection of the adjacency lists of a⃗_x and a_y. Hence, the computation of 2m lists of deltas is already accounted for in Equation <ref>. 




The CBM representation of a binary matrix A is fully characterized by the edges of an MST of the extended distance graph G, and the lists of positive (Δ^+) and negative (Δ^-) deltas associated with each vertex of G. Therefore, the memory required to represent a matrix in CBM format depends on the number of rows in A and the sum of the size of all lists of deltas.



The space required to represent a binary matrix 𝐀∈{0,1}^m × n in Compressed Binary Matrix (CBM) format is proportional to O(m + ∑_i=0^m (|Δ^+_x, r_x| + |Δ^-_x, r_x|)), where r_x represents the index of the row selected to compress row a⃗_x.


    
        Assuming matrix 𝐀 is represented in CBM format. Then the chain of compression required to represent this matrix consists on: (1) a list of edges that represents any MST rooted in vertex 0 of graph G, and (2) a list of positive and negative deltas for each edge of this MST. Since the extended version of G is a fully-connected graph with m+1 vertices, it is known that any MST of G contains m edges. Therefore, the list of edges contains m elements, and there are m lists of deltas, whose size totals ∑_x=0^m (|Δ^+_x,r_x| + |Δ^-_x,r_x|).
    




 §.§ Fast Matrix Multiplication with Compressed Binary Matrix Format


Let w∈ℝ^n be a dense and real-valued vector, and a⃗_x and a⃗_y two distinct rows of a binary matrix 𝐀 as previously defined. 
It follows from Equation <ref> that we can resort to the inner-product a⃗_y ·w⃗ to compute a⃗_x ·w⃗ as

    a⃗_x ·w⃗  = ((a⃗_y ∪Δ^+_x,y )∖Δ^-_x,y)·w⃗ = (a⃗_y ·w⃗) + (Δ^+_x,y·w⃗) - (Δ^-_x,y·w⃗).

This implies that the dot-product a⃗_x ·w⃗ can be calculated in 1+|Δ^+_x,y|+|Δ^-_x,y| scalar operations, if the value of a_y·w⃗ can be reused. Naturally, we can resort to this strategy to design fast matrix-vector multiplication kernels u⃗𝐀v⃗, where 𝐀∈{0,1}^m × n, u∈ℝ^m, and v∈ℝ^n, by computing all dot-products between the rows of 𝐀 and v in an order where: (1) the value of the dot-product a⃗_y ·v is known before calculating a⃗_x ·v, and (2) the value of all dot-products a⃗_x ·v⃗ is calculated with respect to the value a⃗_y ·v⃗ that results in the minimum overall number of scalar operations. By definition, the chain of compression of matrix 𝐀 already represents such an ordering. Therefore, we can accelerate matrix-vector multiplication by traversing the chain of compression of matrix 𝐀 in topological order, and for each edge visited compute u_xa_x ·v as


    u_x  u_y + (Δ^+_x,y·v⃗) - (Δ^-_x,y·v⃗),


where u_y is known to already contain the value a_y ·v. Note that classic sparse-dense matrix-vector multiplication kernels compute u𝐀·v in 𝐧𝐧𝐳(𝐀) scalar operations, where 𝐧𝐧𝐳(𝐀) represents the numbers of nonzero elements in matrix 𝐀. Assuming 𝐀 is represented in CBM format, then the number of deltas required to represent each row a⃗_x of 𝐀 is known to be smaller than, or equal to, the number of nonzero elements in a⃗_x. If the number of deltas required to represent a⃗_x is strictly smaller than 𝐧𝐧𝐳(a⃗_x), then it is clear that the dot-product a⃗_x·v⃗ requires at most 𝐧𝐧𝐳(a⃗_x) scalar operations. On the other hand, if the number of deltas required to represent a⃗_x is the same as the number of nonzero elements in a⃗_x, Equation <ref> suggests that the dot-product a⃗_x ·v⃗ would be computed in 1+𝐧𝐧𝐳(a⃗_x) scalar operations. This scenario can be avoided by engineering the MST algorithm to select an out-going edge of the virtual node 0 in case of draw. Since a⃗_x was not compressed in this scenario, the number of scalar operations required to compute this dot-product is exactly 𝐧𝐧𝐳(a⃗_x).  comment: This next part is confusing. Maybe we could just say that if the nnz of a_x is equal to the deltas then this edge is not considered during the construction of the formatHowever, if the number of deltas required to represent row a⃗_x is the same as the number of nonzero elements in a⃗_x, then this row can be compressed with respect to the null-row a⃗_0. This implies that the dot-product a⃗_x ·v⃗ can be computed in exactly 𝐧𝐧𝐳(a⃗_x), since a⃗_0·v=0, meaning that the cost of this addition can be engineered away. As the number of scalar operations required to compute a_x ·v never surpasses the number of nonzero elements in a⃗_x for x=1,…,m, the following property becomes evident:


    The number of scalar operations required to compute matrix-vector multiplication based on the Compressed Binary Matrix (CBM) format is never greater than those required to compute matrix-vector multiplication based on classic sparse formats.


Additionally, Equation <ref> shows that our matrix-vector multiplication strategy does not require the allocation of additional buffers, since the value of the dot-product required to compute a_x ·v is guaranteed to already be stored in vector u⃗. Therefore, the following property is observed:


    The amount of memory required to compute matrix-vector multiplication based on the Compressed Binary Matrix (CBM) format is proportional to the size of its operands and remains constant during execution time. 


Intuitively, we can resort to the matrix-vector multiplication strategy described above to design fast matrix-matrix multiplication kernels as described in Algorithm <ref>. This algorithm assumes that the left-hand side operand matrix 𝐀∈{0,1}^m× n is represented in CBM format, while matrices 𝐁∈ℝ^n × p and 𝐂∈ℝ^m× p are dense and correspond to the right-hand side operand matrix and to the product matrix, respectively. As it can be observed, Algorithm <ref> computes the matrix-vector product between matrix 𝐀 and each column of matrix 𝐁. Therefore, we can conclude that Properties <ref> and <ref> hold for the matrix-matrix multiplication case.




  
Leveraging High Performance SpMM kernels.  The representation of a binary matrix 𝐀 in CBM format was until now conceptualized as a chain of compression, where each node of this chain is associated with two lists of deltas. As is, a compressed matrix cannot be represented in a sparse format capable of leveraging efficient SpMM kernels. To address this issue, we represent the lists of deltas that characterize our format as a matrix 𝐀'∈ℝ^m × n, which can be represented in any convenient matrix format, and leverage efficient SpMM kernels provided by Intel MKL to compute all dot-products of Algorithm <ref> in a single matrix product 𝐀'𝐁. Once 𝐀'𝐁 is stored in matrix 𝐂, we can finalize the matrix multiplication with CBM, by traversing the chain of compression in topological order, and updating row c⃗_x of 𝐂 as c⃗_x := c⃗_x + c⃗_y for each edge (y,x) that was visited. Naturally, these updates can also leverage efficient AXPY kernels, which are also provided by Intel MKL.



  
Leveraging High Performance SpMM kernels.  The representation of a binary matrix 𝐀 in CBM format was until now conceptualized as a chain of compression, where each node of this chain is associated with two lists of deltas. As is, a compressed matrix cannot be represented in a sparse format capable of leveraging efficient SpMM kernels. To address this issue, we represent the lists of deltas that characterize our format as a matrix 𝐀'∈ℝ^m × n, which can be represented in any convenient matrix format, and leverage efficient SpMM kernels provided by Intel MKL to compute all dot-products of Algorithm <ref> in a single matrix product 𝐀'𝐁. Once 𝐀'𝐁 is stored in matrix 𝐂, we can finalize the matrix multiplication with CBM, by traversing the chain of compression in topological order, and accordingly updating the value of row c⃗_x of 𝐂 with the current value of c⃗_y. Naturally, these updates can also leverage efficient AXPY kernels, which are also provided by Intel MKL.



  
Multi-threading parallelism. As suggested in the previous paragraph, matrix multiplication with the CBM format can be divided into two stages. The first one computes the product 𝐂:=𝐀'𝐁 which is embarrassingly parallel, and Intel MKL already provides efficient multi-threaded and vectorized implementations of this operation. The second stage involves updating the rows of matrix 𝐂 with respect to the chain of compression. This stage presents data-dependencies, since the final value of row c⃗_x of 𝐂 can only be calculated once c⃗_y is known. There are however no dependencies between different branches of the chain of compression. Therefore, we can parallelize this stage 
by concurrently updating the rows of matrix 𝐂 that are found in different branches of the virtual node 0.



  
Extending SpMM with CBM to normalized adjacency matrices. Let 𝐀̂∈ℝ^n × n represent the normalized adjacency matrix of an unweighted graph, where 𝐀̂ = 𝐃^-1/2 (𝐀 + 𝐈) 𝐃^-1/2. As is, we cannot resort to the CBM format to represent 𝐀̂ since this matrix is not binary. However, fast matrix multiplication is still possible. Note that 𝐃^-1/2 is a diagonal matrix, and therefore (𝐀 + 𝐈) 𝐃^-1/2 corresponds to a column-scaled matrix, i.e., the elements in a column of the matrix are either 0 or have a constant value that is unique to this column. We can represent this matrix in CBM format, by simply multiplying the corresponding matrix of deltas (𝐀' + 𝐈) by 𝐃^-1/2. At this point, we can efficiently compute 𝐂:=((𝐀 + 𝐈)𝐃^-1/2)𝐁 as described previously. 
Finally, note that 𝐃^-1/2(𝐀 + 𝐈)𝐃^-1/2𝐁), simply scales the rows of matrix (𝐀 + 𝐈)𝐃^-1/2𝐁. The cost of scaling the rows of this matrix can be hidden by fusing it with the update step of the product matrix 𝐂. 



  
Speeding up SpMM with Edge Pruning. Note that not all compression opportunities contribute to faster matrix multiplication kernels. The overheads associated with differential compression, such as traversing the chain of compression and updating the result matrix, might overcome potential performance gains if the number of scalar operations saved are not above a certain threshold. To address this issue, we can prune all edges of the distance graph of 𝐀 where the number of scalar operations saved does not meet a user-defined threshold α∈ℕ. For each edge (y,x) in the extended distance graph of 𝐀, prune this edge if its weight is greater than the number of nonzero elements in a⃗_x minus α. Naturally, if we prune the edges of the extended distance graph of 𝐀 in this manner, it is possible for a single edge direction to be pruned, while the opposite direction remains in the graph. Therefore, the extended distance graph of 𝐀 is now directed, and a suitable chain of compression corresponds to a Minimum Cost Arborescence (MCA) rooted in the virtual node 0. Note that our compression algorithm remains correct, since the extended distance graph of 𝐀 contains an out-going edge from the virtual node 0 to all other nodes; Furthermore, the time required to build the CBM format, as shown in Lemma <ref>, remains unchanged since finding an MCA or an MST present the same time complexity  <cit.>.



§ EXPERIMENTAL EVALUATION



  
Experimental Setting. The experiments found in this section were run on an Intel Xeon Gold 6130 (Skylake) CPU with 16 physical cores and 2.1 GHz fixed clock frequency. This machine runs on CentOS Linux 7 (version 3.10.0) operating system. The SpMM kernels tested in this section were implemented in C++, and rely on Intel MKL 24.0 sparse CSR format and corresponding matrix multiplication kernels. The C++ code developed in this work is then called from Python 3.11, via PyTorch C++ Extensions to reliably emulate common use-cases. 
Parallel experiments (with 16 cores) were implemented with OpenMP 4.5, and the threads were pinned to physical cores with  environment variable.




    
    
    
    
    






 §.§ Evaluation Metrics


The CBM format was evaluated with respect to quality of compression achieved (compression ratio), and the time required (runtime reduction) to compute SpMM and to infer a 2-layer GCN, when the adjacency
matrix of the graphs are represented in our format. The compression
ratio is defined as the ratio between the memory required to represent a matrix in CSR format and 
CBM format. In our implementation, the CBM format is composed by the corresponding matrix of deltas 𝐀' and a tree representing the chain of compression, both stored in CSR format. In the context of sparse-dense matrix multiplication (SpMM), the runtime reduction is measured by comparing the average time, out of
50 runs, it takes to perform matrix multiplication with a randomly generated
dense matrix with 500 columns using the CSR format, to the time taken to compute
the same matrix multiplication with the CBM format. The formula used to capture
this metric is T_CSR - T_CBM/T_CSR×
100%, where T_CSR is the time required to carry out sparse-dense matrix multiplication with CSR by the state-of-the-art SpMM implementation offered by Intel MKL, and T_CBM is the time taken to compute the same product with CBM. We used same formula and number of runs in the context of GCN inference, however, T_CSR and T_CBM correspond to the time required by the inference stage of this network by resorting to SpMM kernels based on CSR and CBM, respectively. It is important to note that we did not consider the SpMM kernels that are native to PyTorch, because these kernels were substantially slower than the ones implemented in Intel MKL.  



 §.§ Datasets


To demonstrate the advantages of the CBM format in the context of SpMM and GCN inference, we selected six real-world graphs of varying size and average in-degree as depicted in Table <ref>. The selected graphs depict relationships between authors and/or academic papers, where nodes tend to share many common neighbors. This property suggests that the adjacency matrices of these graphs are good candidates to be represented in CBM format. 

In co-paper graphs, each node represents a paper. An undirected edge is placed between two nodes if the corresponding papers share at least one common author. They depict the interconnection and collaborative patterns between various academic publications.

Co-author graphs represent scientific collaborations between authors of academic papers. Here, nodes correspond to authors, and an undirected edge is placed between two nodes if authors of the nodes have co-authored a paper together. If a paper is authored collaboratively by a group of authors, it results in a fully connected subgraph, or clique, encompassing those grouped nodes.

Citation graphs are directed graphs where each node represents an academic paper. Directed edges in these graphs illustrate citations, with an edge pointing from the citing paper to the cited paper. These graphs highlight the directional flow of information and the influence of one paper upon another within the academic community.







 §.§ Sparse-Dense Matrix Multiplication (SpMM) Evaluation
 

Finding the best α is key to improve the performance of matrix multiplication with CBM. Adjusting this parameter not only reduces overhead associated with
traversing the compression chain, but also exposes more parallelism opportunities as it increases the out-degree of the virtual node. Given the importance of α we first consider the case where α = 0 and our edge pruning technique was not applied, and then we show how fine-tuning α improves the quality of matrix multiplication with CBM.   







impacts the The behaviors of α for each datasets when setting the α parameter to different values, impacting both compression effectiveness and computational performance. This paragraph is very confusing. Please be direct, no point on trying to be fancy



  
α equal to zero. Representing Cora (Fig. <ref>) and PubMed (Fig. <ref>) datasets in CBM resulted in minimal, and even negative, compression
gains with respect to CSR. This is likely caused by the small average in-degree of these graphs, which suggests that the compression opportunities found in these graphs do not offset the memory overhead required to represent the corresponding chains of compression. As expected, the poor compression rate of these datasets led to no speedup in the context of matrix multiplication. 
On the other hand, compressing ca-AstroPh (Fig. <ref>) and ca-HepPH (Fig. <ref>) with our format respectively increased the compression ratios of both datasets to 2.6× and 1.7×, subsequently accelerating matrix product for both datasets. Our matrix multiplication strategy achieved a runtime reduction of 26% and 8% for ca-AstroPh in sequential and parallel environments, respectively. For ca-HepPh the same multiplication kernel presented a runtime reduction of 40% and 18% also in sequential and parallel environments. 
The datasets that benefited the most from our format were the coPapersCiteseer (Fig. <ref>) and coPapersDBLP (Fig. <ref>), achieving impressive compression ratios of 9.8× and 5.9×, leading to substantial improvements in computational performance. Our matrix multiplication kernel with coPapersCiteseer achieved a runtime reduction of 71% and 77% in sequential and parallel environments, while the same kernel with coPapersDBLP presented a runtime reduction of 59% and 61% in the same experimental settings. These results highlight the potential of the CBM format to efficiently compress and multiply unweighted graphs that present natural communities and an high average in-degree.






  
α greater than zero. Setting α to a value greater than 8 reduced the overhead associated with traversing the chain of compression of Cora (Fig. <ref>) and PubMed (Fig. <ref>), making the performance decay observed for α = 0 negligible in sequential and parallel settings, even when no compression gains are observed. 
Increasing the value of α is also beneficial for sequential matrix multiplication with ca-AstroPh (Fig. <ref>) and ca-HepPh (Fig. <ref>), improving the respective runtime reduction from 26% up to 28% (α=2) and from 40% up to 44% (α=4). Nevertheless, it is important to note that our compression algorithm will start to ignore good compression opportunities once α is large enough, worsening the performance of our matrix multiplication strategy. This effect is evident in both datasets for α greater than 16, where the compression ratio decreases alongside with the runtime reduction. 
Experiments with both co-author datasets show that adjusting α is even more important in the parallel case, increasing the runtime reduction of our multiplication strategy with ca-AstroPh from 8% up to 11% (α=8), and with ca-HepPh from 18% up to 26% (α=64). Furthermore, these experiments confirm that the degree of parallelism of our multiplication strategy increases concurrently with α. This effect is easily observed for ca-HepPh (Fig.<ref>), where the performance of our matrix multiplication kernel sharply declines for α greater than 2, followed by a steep increase in performance when α equals to 32 or 64 (even though our compression algorithm is already ignoring good compression opportunities at this point). 
Finding the best α is not as relevant for coPapersCiteseer (Fig.<ref>)  and coPapersDBLP (Fig.<ref>), as most compression opportunities save more than 8 scalar operations in the context of our matrix multiplication strategy. This observation is verified, since the compression ratio of these datasets shows little to no decrease for α smaller than 8. This effect is most likely due to the high average in-degree of both datasets. Still, adjusting α is required to obtain the best runtime reduction for both coPapers datasets. In our experiments our matrix multiplication kernel with coPapersCiteseer achieved a peak runtime reduction of 71% (α=2) and 79% (α=32) in sequential and parallel environments, while the same kernel with coPapersDBLP peaks at 59% (α=4) and 63% (α=16) also in sequential and parallel environments.


increasing α also tends to increase the degree of parallelism of our matrix multiplication strategy. This effect can be observed    

 Experiments in a parallel environment Adjusting α also increases  ca-Astro-Ph and ca-HepPh  we can observe that increases the degree of parallelism 


For the same datasets, adjusting α in the context of parallel matrix multiplication resulted in a runtime reduction up to 11% (α=8) and 26% (α=64). 



Increasing the value of α is also beneficial for matrix multiplication with ca-AstroPh (Fig. <ref>) and ca-HepPh (Fig. <ref>), respectively improving the sequential runtime reduction of these datasets up to 28% (α=2) and 44% (α=4). For the same datasets, adjusting α in the context of parallel matrix multiplication resulted in a runtime reduction up to 11% (α=8) and 26% (α=64). 

While increasing α tends to improve the performance of our matrix multiplication strategy, it is important to note that our compression algorithm will start to ignore good compression opportunities once α becomes large enough, thus worsening the performance of our matrix multiplication strategy. This effect is evident in the context of sequential matrix multiplication for α greater than 16 in both datasets.

Fine-tuning α in a parallel environment is more complex. As it can be observed in Figure <ref>, increasing α represents a trade-off between good compression ratio and improving the degree of parallelism of our solution. This effect is easily observed for matrix products involving ca-HepPh (Fig.<ref>), where the performance of our matrix multiplication kernel sharply declines for α greater than 2, followed by a steep increase in performance when α equals to 32 or 64. This increase in performance is the result of         


 




For ca-AstroPh, as α increases, the sequential runtime initially
improves, peaking at 28% reduction for α = 2 before dropping
significantly to just 4% for α = 64. Similarly, the parallel runtime
initially increases to 11% for α = 8 and then declines to 2% for
α = 64. The compression ratio decreases steadily, reaching 1.02×
for α = 64. For ca-HepPH, the sequential runtime improvement peaks at
44% for α = 4, then declines to 29% for α = 64, while the
parallel runtime peaks at 19% for α = 1  and eventually increases again
to 26% for α = 64. The compression ratio gradually falls to
1.56× as α increases.

For coPapersDBLP, as α increases, the sequential runtime peaks slightly
at 59% for α = 4 before reducing further to 42% for α = 64, and
the parallel runtime increases marginally before dropping to 59% for α =
64. The compression ratio declines to 2.2×.

In coPapersCiteseer, the sequential runtime improves marginally before
decreasing to 59% for α = 64, and parallel runtime slightly improves,
peaking at 79% for α = 32 before slightly adjusting to 77% for α
= 64, with the compression ratio falling to 3.4×.

While the compression rate decreases as α increases, we observe better
SpMM performance overall. The impact of a better α becomes more
critical in the parallel case because it enhances the degree of parallelism of
the SpMM kernel. This is evident as in some cases the performance decreases as
α increases, and as we further increase the value of alpha, the
performance spikes. For example, ca-HepPh shows a decrease from α=8 to
α=16, but spikes from α=16 to α=32. If α is large
enough (around α=8), we can minimize the performance decay of
challenging datasets such as Cora or PubMed.






 §.§ Integrating CBM with GCN


To assess the impact of CBM in GNNs, we considered the runtime reduction of an inference stage of a 2-layer GCN with 500 features where the normalized adjacency matrix 
of each dataset is represent in our format. Matrix products involving the normalized adjacency matrix 
were carried out with our extended matrix multiplication kernel, as described in Section <ref>. As baseline for our experiments, we selected the same 2-layer GCN where the normalized adjacency matrix 
is represented in CSR and any products involving this matrix are carried out by the SpMM implementation found in Intel MKL.
These experiments are illustrated in Figure <ref>. To keep the discussion concise, we only considered the values of α that led to the best matrix multiplication performance for each dataset analyzed in Section <ref>.

Representing both Cora and PubMed datasets in CBM format increased the inference time of the corresponding GCNs compared to the baseline. This behavior is expected because the only steps of the inference that benefit from the CBM format are the products involving the normalized adjacency matrix. As previously shown, compressing these datasets using our format did not accelerate the corresponding matrix products.
Experiments with both co-author datasets demonstrated that our format can reduce the inference time for GCN models. For ca-AstroPh, our format reduced the inference time of the network by 17% in sequential environments and 3% in parallel environments. Similarly, for ca-HepPh, our format achieved a runtime reduction of 22% in sequential environments and 11% in parallel environments.
Representing both co-Paper datasets in our format resulted in the highest runtime reductions during GCN inference. For coPapersCiteseer, our format achieved an average runtime reduction of 48% in sequential environments and 66% in parallel environments. For coPapersDBLP, our format achieved a runtime reduction of 43% in sequential environments and 52% in parallel environments.


Significant reductions in inference time for
coPapersDBLP and coPapersCiteseer by 43% and 52% in sequential mode, and 48%
and 66% in parallel mode, respectively,

for each datasetThe baseline of our experiments is the same 2-layer GCN, but the normalized adjacency matrix is represented in CSR. Products involving this matrix were carried out by resorting to the SpMM implementation offered by Intel MKL. For the sake of briefness, the experiments of Figure <ref> only consider the α va for each dataset, as seen in Section <ref>) that resulted in the best runtime reduction for both sequential and parallel execution. 

As part of the inference stage of a GCN, we integrated CBM SpMM
for normalized adjacency matrices to increase efficiency.
Figure <ref> illustrates the average runtime reduction for the inference
time. We set the baseline against the CSR SpMM provided by Intel MKL for each
dataset. Runtime was measured by the average of 50 epochs with two GCN
layers. Based on this comparison, it can be seen that CBM can enhance
computations within GCN architectures.


Cora and PubMed datasets show an increase in inference time by 6% and 7% in
sequential mode, and 19% and 8% in parallel mode, respectively, indicating
scenarios where CBM format might not be optimal. In contrast, ca-AstroPh and
ca-HepPh datasets benefit from reductions in inference time by 17% and 3% in
sequential mode, and 22% and 11% in parallel mode, respectively, suggesting
moderate compatibility with CBM.  Significant reductions in inference time for
coPapersDBLP and coPapersCiteseer by 43% and 52% in sequential mode, and 48%
and 66% in parallel mode, respectively, demonstrate superior performance of CBM
SpMM for datasets favorable to block compression.



§ FINAL REMARKS


In this work we proposed the Compressed Binary Matrix (CBM) format which simultaneously reduces the memory footprint of unweighted graphs and binary matrices, and enables the implementation of new matrix multiplication kernels that might be significantly faster than the current state-of-the-art.
Experimental evaluation results shown significant speedups, both in sequential and parallel environments, up to 5×.
We obtained also significant performance improvements in the context of the GCN inference stage by integrating the CBM format in a deep learning framework, namely PyTorch, observing speed ups of 3×.
Although we did not discuss the CBM format construction time, we observe that it can be built in a reasonable amount of time for a dataset provider. In our experiments it took us less than 16 seconds to convert the largest dataset into our format in a sequential CPU environment.
It is important to stress that the effectiveness of our format depends on the specific dataset, as discussed in the experimental evaluation. While we suspect that graphs with a high average degree and a tendency to form communities are good candidates, the best way to determine if a graph is suitable is to examine the compression ratio achieved by our format for a reasonable value of α. Finally we highlight that our format is future proof, since future optimizations to high-performance SpMM kernels will also accelerate matrix multiplication with the CBM format.

Future work concerns integrating and evaluating the CBM format in the context of different GNNs architectures, and also targeting the training stage of this networks. Additionally, we intend to implement and evaluate our format and corresponding multiplication kernels in GPU architectures. 



§ ACKNOWLEDGMENTS

This work has been supported by the Innovation Study CBM4scale, funded by the Inno4scale project, which is funded by the European High-Performance Computing Joint Undertaking (JU) under Grant Agreement No 101118139. The JU receives support from the European Union's Horizon Europe Programme. 



§ AUTHOR CONTRIBUTIONS


JNFA devised the main conceptual of the Compressed Binary Matrix (CBM), including the various matrix multiplication algorithms and their optimizations. JNFA also implemented the CBM format and the related multiplication kernels in C++, including the integration of SpMM and AXPY (from Intel MKL) into the matrix multiplications kernels based on the CBM format, as proposed in Section <ref>. Additionally, JNFA implemented an interface that enables calling the previous C++ routines from Python.
SM designed the Python benchmarks to compare CBM and CSR in matrix multiplication and compression quality, as described in Section <ref>. SM also integrated the CBM and CSR-based matrix multiplication kernels into the Message Passing Layer (MPL) to evaluate the impact of the CBM format during the inference stage of a 2-layer GCN model implemented in PyTorch, as discussed in Section <ref>. SB, APF, WNG, and LMSR supervised the project. JNFA wrote the bulk of this draft. All authors provided critical feedback, participated in discussions, contributed to the interpretation of the results, and approved the final manuscript.


unsrt  
  


