#  SAGE: A Storage-Based Approach for Scalable and Efficient Sparse
Generalized Matrix-Matrix Multiplication

Myung-Hwan Jang  Department of Computer Science
Hanyang UniversitySeoulRepublic of Korea<EMAIL_ADDRESS>0000-0003-4419-5148 ,  Yunyong Ko11footnotemark: 1 Department of Computer
Science
Hanyang UniversitySeoulRepublic of Korea<EMAIL_ADDRESS>0000-0003-1283-4697 ,  Hyuck-Moo Gwon  Department of Computer Science
Hanyang UniversitySeoulRepublic of Korea<EMAIL_ADDRESS>0000-0002-8765-1863 ,  Ikhyeon Jo  Department of Computer Science
Hanyang UniversitySeoulRepublic of Korea<EMAIL_ADDRESS>0000-0003-2245-6220 ,  Yongjun Park22footnotemark: 2 Department of Computer
Science
Yonsei UniversitySeoulRepublic of Korea<EMAIL_ADDRESS>0000-0003-3725-0380 and  Sang-Wook Kim  Department of Computer Science
Hanyang UniversitySeoulRepublic of Korea<EMAIL_ADDRESS>0000-0002-6345-9084

(2023)

###### Abstract.

Sparse generalized matrix-matrix multiplication (SpGEMM) is a fundamental
operation for real-world network analysis. With the increasing size of real-
world networks, the single-machine-based SpGEMM approach cannot perform SpGEMM
on large-scale networks, exceeding the size of main memory (i.e., not
scalable). Although the distributed-system-based approach could handle large-
scale SpGEMM based on multiple machines, it suffers from severe inter-machine
communication overhead to aggregate results of multiple machines (i.e., not
efficient). To address this dilemma, in this paper, we propose a novel
storage-based SpGEMM approach (SAGE) that stores given networks in storage
(e.g., SSD) and loads only the necessary parts of the networks into main
memory when they are required for processing via a 3-layer architecture.
Furthermore, we point out three challenges that could degrade the overall
performance of SAGE and propose three effective strategies to address them:
(1) block-based workload allocation for balancing workloads across threads,
(2) in-memory partial aggregation for reducing the amount of unnecessarily
generated storage-memory I/Os, and (3) distribution-aware memory allocation
for preventing unexpected buffer overflows in main memory. Via extensive
evaluation, we verify the superiority of SAGE over existing SpGEMM methods in
terms of scalability and efficiency.

sparse matrix multiplication, real-world graphs, network analysis

††journalyear: 2023††copyright: acmlicensed††conference: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management;
October 21–25, 2023; Birmingham, United Kingdom††booktitle: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management
(CIKM ’23), October 21–25, 2023, Birmingham, United Kingdom††price:
15.00††doi: 10.1145/3583780.3615044††isbn: 979-8-4007-0124-5/23/10††ccs:
Information systems Data management systems

##  1\. Introduction

Graphs are widely used to model real-world networks, where a node represents
an object and an edge does the pair-wise relationship between two objects. To
analyze real-world networks and discover useful knowledge, many graph
algorithms have been studied (Sedgewick and Wayne, 2011; Liu et al., 2016;
Page et al., 1999; Kleinberg, 1999; Suzuki et al., 2003; Karypis and Kumar,
1998; Becchetti et al., 2008; He et al., 2010; Kepner and Gilbert, 2011;
Kepner et al., 2016; Kang et al., 2020; Hamedani and Kim, 2017; Ko et al.,
2021b), where a graph is generally represented as a dense matrix (i.e., 2-D
array), and each element $e_{i,j}=1$ if an edge exists between nodes $i$ and
$j$, and $e_{i,j}=0$ otherwise. Generalized matrix-matrix multiplication
(GEMM), one of the key operations in these algorithms, plays a fundamental
role to diffuse the node information via network connectivity. Recently, with
the growing size of real-world networks, it has become more crucial to
efficiently perform GEMM on large-scale graphs (Gu et al., 2017; Zheng et al.,
2016; Azad et al., 2016).

$10^{1}$$10^{2}$$10^{3}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$$10^{5}$Node
degreeNumber of nodes (a) Degree distribution.

(b) Matrix density ($\approx 0.15\%$).

Figure 1. (a) Power-law degree distribution and (b) visualized matrix density
of real-world graphs.

Performing GEMM on real-world graphs represented as dense matrices, however,
often requires an unnecessarily large amount of space and computational cost
(Nagasaka et al., 2018; Pal et al., 2018). This is because they tend to follow
a power-law degree distribution (Leskovec et al., 2007) – i.e., a majority of
nodes have only a few edges, while a small number of nodes have a large number
of edges – so that the number of non-zero elements (existing edges) is much
smaller than that of zero elements (non-existing edges) in the matrix for a
real-world graph. Figure 1(a) shows the power-law degree distribution of
Wikipedia (Jo et al., 2019) and Figure 1(b) does the matrix density of Wiki-
Vote (Leskovec et al., 2010b, a), where each point represents a non-zero
element, respectively.

To address this issue, sparse GEMM (SpGEMM) methods have been widely studied
(Akbudak and Aykanat, 2017; Demirci and Aykanat, 2020; Azad et al., 2016;
Nagasaka et al., 2017; Deveci et al., 2017; Gremse et al., 2018; Anderson et
al., 2016; Pal et al., 2018; Zhang et al., 2020; Davis, 2019), where a graph
is represented as a sparse matrix containing only existing edges. On one hand,
single-machine-based methods (Akbudak and Aykanat, 2017; Wang et al., 2014;
Xie et al., 2019; Nagasaka et al., 2018) aim to perform SpGEMM efficiently on
a single machine by considering the characteristics of real-world graphs,
where they assume that the sparse matrix can fit in main memory. Thus, they
are unable to process SpGEMM on large-scale graphs exceeding the size of the
main memory (i.e., not scalable). On the other hand, distributed-system-based
methods (Buluç and Gilbert, 2012; Azad et al., 2016; Selvitopi et al., 2019;
Demirci and Aykanat, 2020) are able to perform SpGEMM on large-scale graphs,
exceeding the size of the main memory, by using multiple machines
concurrently. These methods, however, require a substantial amount of inter-
machine communication overhead to aggregate the results from multiple
machines, which is often more than $50\%$ of the entire process (Azad et al.,
2016; Ko et al., 2021a), and a lot of costs and efforts to maintain complex
distributed systems (i.e., not efficient).

To tackle these two limitations together, in this paper, we propose a novel
storage-based approach for large-scale SpGEMM, named Storage-bAsed approach
for SpGEMM (SAGE). SAGE stores the entire graphs in external storage of a
single machine and loads only some parts of the graphs when necessary into the
main memory to process them. For efficient data transfers between storage and
main memory, we design SAGE with a 3-layer architecture, i.e., (1) storage,
(2) in-memory, and (3) operation layers, where each layer closely interacts
with other layers. Thus, SAGE (1) is able to process large-scale graphs,
exceeding the size of the main memory by using sufficient capacity of external
storage (i.e., scalable) and (2) requires only intra-machine communication
overhead (storage-memory I/Os), much smaller than the inter-machine
communication overhead (i.e., efficient).

To improve the performance of SAGE, it is crucial to handle storage-memory
I/Os efficiently, which are inevitably generated due to the limited size of
main memory on a single machine. To this end, first, we point out three
important challenges that could lead to serious performance degradation in
SAGE: (C1) Workload balancing, i.e., How to distribute workloads of SpGEMM
evenly across multiple threads?, (C2) Intermediate handling, i.e., How to
handle a large amount of the storage-memory I/Os generated by the intermediate
results of SpGEMM?, and (C3) Memory management, i.e., How to allocate the main
memory space to three buffers to reduce storage-memory I/Os? (for two input
and one output matrices). Then, we propose three effective strategies to
address these challenges: (1) block-based workload allocation to evenly
distribute workloads of SpGEMM into multiple threads for C1, (2) in-memory
partial aggregation to reduce the amount of unnecessarily generated storage-
memory I/Os for C2, and (3) distribution-aware memory allocation to adjust the
proportions of the three buffers, based on the characteristics of real-world
graphs for C3.

We note that although there exists several existing storage-based methods
(Kyrola et al., 2012; Roy et al., 2013; Han et al., 2013; Zheng et al., 2015;
Jo et al., 2019), they focus only on the multiplication between sparse matrix
(SpM) and dense vector (DV), i.e., SpMV but not SpGEMM that we focus on. More
specifically, the existing storage-based SpMV methods handle only “one” large
sparse matrix and small dense vectors ($SpM\times DV=DV$), where the dense
vectors are much smaller than a sparse matrix ($SpM>>DV$). On the other hand,
SAGE handles “three” large sparse matrices ($SpM\times SpM=SpM$), which
implies that it is more challenging to efficiently process storage-memory I/Os
in the storage-based SpGEMM than the storage-based SpMV. To the best of our
knowledge, this is the first work to successfully perform large-scale SpGEMM
by using external storage on a single machine. We believe that SAGE could be a
practical alternative to researchers and practitioners who aim to analyze
large-scale real-world graphs.

The main contributions of this work are summarized as follows:

  * •

New Approach: We propose a new storage-based approach for large-scale SpGEMM,
SAGE to address the limitations of existing single-machine-based (scalability)
and distributed-system-based (efficiency) approaches simultaneously.

  * •

Effective Strategies: We point out three performance-critical challenges and
propose effective strategies to tackle them: (1) block-based workload
allocation, (2) in-memory partial aggregation, and (3) distribution-aware
memory allocation.

  * •

Extensive Evaluation: We conduct extensive evaluation, which demonstrates that
(1) SAGE is able to perform SpGEMM on large-scale graphs, surpassing the
limits of the existing single-machine-based methods, (2) SAGE achieves high
SpGEMM performance comparable to or better than the distributed-system-based
methods, and (3) each of the proposed strategies of SAGE is effective in
improving the performance of SAGE.

##  2\. Related Work

In this section, we review existing SpGEMM methods. Table 1 compares SAGE with
existing methods in terms of the three performance-critical challenges: (C1)
workload balancing (WB), (C2) intermediate handling (IH), and (C3) memory
management (MM).

Table 1. Comparison of existing SpGEMM methods with SAGE based on the three performance-critical challenges Method | (C1) WB | (C2) IH | (C3) MM
---|---|---|---
Intel MKL (Wang et al., 2014) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  |
KNL-SpGEMM (Nagasaka et al., 2018) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  |
IA-SpGEMM (Xie et al., 2019) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  |
SpSUMMA (Buluç and Gilbert, 2012) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  |
Graphulo (Hutchison et al., 2015) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |
Split-3D (Azad et al., 2016) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  |
gRRp (Demirci and Aykanat, 2020) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |
SAGE (proposed) |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓ |  TextRenderingMode=FillStroke, LineWidth=.5pt, ✓

Single-machine-based approach. A single-machine-based approach (Wang et al.,
2014; Xie et al., 2019; Akbudak and Aykanat, 2017; Nagasaka et al., 2018) aims
to improve the performance of SpGEMM, based on the properties of real-world
networks on a single machine (e.g., node degree distribution). Intel MKL (Wang
et al., 2014) is an open source library supporting a series of math functions
such as SpGEMM and SpMV, which are optimized for Intel CPU architecture.
Nagasaka et al. (2018) proposed memory management to keep output results and
thread scheduling strategies to resolve the bottlenecks of the entire SpGEMM
process. Xie et al. (2019) proposed an input-aware auto-tuning framework for
SpGEMM (IA-SpGEMM) to reduce the overhead of memory access and sparse
accumulation in SpGEMM. These single-machine-based methods, however, assume
that the entire input matrices can be loaded in the main memory (i.e., in-
memory based SpGEMM) because the input matrices have only existing edges in
SpGEMM (i.e., sparse matrix). Therefore, they are not able to process SpGEMM
on large-scale graphs exceeding the size of main memory in a single machine.

There are some works using external storage to process large-scale graphs on a
single machine (Kyrola et al., 2012; Han et al., 2013; Zheng et al., 2015,
2016; Jo et al., 2019). However, they focus on improving sparse matrix and
dense vector multiplication (SpMV) rather than SpGEMM that our work focuses
on. Note that we can obtain the exactly same results as SpGEMM by iteratively
performing SpMVs in concept; however, this way is not suitable for processing
SpGEMM since the operations for dense vectors can lead to serious performance
degradation due to a large amount of memory usage and computation overhead
(Pal et al., 2018; Buluç and Gilbert, 2012; Azad et al., 2016).

Distributed-system-based approach. A distributed-system-based approach (Buluç
and Gilbert, 2012; Hutchison et al., 2015; Azad et al., 2016; Selvitopi et
al., 2019; Demirci and Aykanat, 2020) aims to process SpGEMM on large-scale
graphs that cannot be loaded on the main memory of a single machine, via
multiple machines in a distributed cluster. This approach splits and stores
input matrices into multiple machines, and then processes them in parallel.
Buluç and Gilbert (2012) proposed SpSUMMA that parallelizes SpGEMM by
partitioning input matrices into 2-D grids and distributing them to multiple
processors. Hutchison et al. (2015) proposed a distributed SpGEMM method,
Graphulo that considers the locality of data stored in a distributed database
to reduce the inter-machine communication overhead. Azad et al. (2016)
proposed Split-3D that splits input/output matrices into 3-D grids for
parallelization of computations and communications. Demirci and Aykanat (2020)
proposed gRRp, a bipartite graph-based partitioning method for balancing
workloads among multiple machines.

Though the distributed-system-based approach is able to perform SpGEMM on
large-scale graphs, a substantial amount of inter-machine communication
overhead is required to aggregate the results from multiple machines, which is
non-trivial (more than $50\%$) (Azad et al., 2016; Ko et al., 2021a). It also
requires a lot of costs and efforts to maintain the high-performance
infrastructure of distributed systems such as fault tolerance and graph
partitioning for distributed machines.

##  3\. SAGE: Proposed framework

Table 2. Notations and their descriptions Notation | Description
---|---
$M^{1}_{in},M^{2}_{in}$ | input matrices of SpGEMM
$M_{out}$ |  output matrix ($M^{1}_{in}\times M^{2}_{in}$)
$M(i,;)$ |  $i$-th row of matrix $M$
$B^{1}_{in},B^{2}_{in}$ | buffers for input matrices in the main memory
$B_{out}$ | buffer for output matrix in the main memory
$T^{1}_{buf},T^{2}_{buf}$ |  buffer index tables for input matrices $M^{1}_{in},M^{2}_{in}$
$T^{1}_{obj},T^{2}_{obj}$ |  object index tables for input matrices $M^{1}_{in},M^{2}_{in}$

We present a novel storage-based SpGEMM method, named Storage-bAsed approach
for SpGEMM (SAGE). First, we describe the notations and the problem statement
(Section 3.1). Then, we present the architecture and core algorithms of SAGE
(Section 3.2), and three novel strategies to address performance-critical
challenges (Section 3.3). Finally, we analyze the complexity of SAGE (Section
3.4).

###  3.1. Notations and Problem Statement

####  3.1.1. Notations

Table 2 shows the notations used in this paper. We denote two input sparse
matrices for SpGEMM as $M^{1}_{in}$ and $M^{2}_{in}$, and the output matrix as
$M_{out}$ (i.e., $M^{1}_{in}\times M^{2}_{in}$). In a matrix $M$, $M(i,;)$ is
an $i$-th row, and $M(i,j)$ is a $j$-th element in the $i$-th row, where the
length of each row ($|M(i,;)|$) can vary across rows because $M$ is a sparse
matrix. For clarity, we denote buffers and index tables in the main memory as
$B$ and $T$, respectively. $B^{1}_{in}$ and $B^{2}_{in}$ are the input buffers
to hold necessary parts of the input matrices, and $B_{out}$ is the output
buffer for the intermediate results.

####  3.1.2. Problem statement

Given two input sparse matrices of $M^{1}_{in}$ and $M^{2}_{in}$, SAGE aims to
produce the output sparse matrix $M_{out}$ by performing SpGEMM between the
two sparse matrices. When performing SpGEMM, three different types of products
can be considered: (1) inner product (row-by-column product), (2) outer
product (column-by-row product), and (3) row-wise product (row-by-row
product). Although these three products generate exactly the same result,
$M_{out}$, they require different computation overhead and memory space.
Specifically, the row-wise product requires not only much smaller computation
overhead (e.g., index-matching) than the inner product, but also less memory
usage than the outer products (Srivastava et al., 2020; Baek et al., 2021).
Thus, we adopt the row-wise product for performing SpGEMM in SAGE, following
existing algorithms (Akbudak and Aykanat, 2017; Wang et al., 2014; Xie et al.,
2019; Nagasaka et al., 2018). In the row-wise product, the $i$-th row of the
output matrix is computed by Eq. 1,

(1) $M_{out}(i,;)=\sum_{j=1}^{n}M^{1}_{in}(i,j)\times M^{2}_{in}(j,;),$

where $n$ indicates the number of rows in the input matrices. Figure 2
illustrates the process of the row-wise product computing the $i$-th row of
the output matrix $M_{out}$ in performing SpGEMM. First, every $j$-th element
of the $i$-th row in the former input matrix, $M^{1}_{in}(i,j)$, is multiplied
by its corresponding $j$-th row in the latter input matrix, $M^{2}_{in}(j,;)$,
(highlighted in the same color in Figure 2). Then, the intermediate results
are aggregated to compute the final result of the $i$-th row of the output
matrix, $M_{out}(i,;)$. This process is repeated until all rows of the output
matrix are computed.

Figure 2. The process of row-wise product. Figure 3. Overview of SAGE with the
3-layer architecture (storage, in-memory, and operation layers).

###  3.2. Architecture and Algorithm

SAGE stores the entire graphs in storage (e.g., SSD), loads only the necessary
parts of the graphs into the main memory, and processes them in parallel by
using multiple threads. Thus, storage-memory I/Os are inevitably generated
during the process of storage-based SpGEMM, which could become a bottleneck of
the overall SpGEMM process. To handle the inevitable overhead efficiently, in
this section, we propose a novel 3-layer architecture for efficient data
transfers between storage and main memory (Section 3.2.1), and describe how
SAGE processes SpGEMM based on the architecture (Section 3.2.2).

####  3.2.1. Architecture

As illustrated in Figure 3, SAGE consists of three layers: (1) storage, (2)
in-memory, and (3) operation layers, where each layer closely interacts with
other layers to efficiently process storage-memory I/Os.

Storage layer. This layer is in charge of managing the data stored in external
storage: two input matrices ($M^{1}_{in}$ and $M^{2}_{in}$), an output matrix
($M_{out}$), and intermediate results ($M_{out}(i,;)$). This storage layer
handles read/write requests from the in-memory layer. We define a fixed-size
unit standard I/O as a block and store the input matrices in a set of blocks.
Each block usually contains multiple objects, each of which stores each row of
the input matrices (i.e., a node and its related edges). In a real-world graph
following the power-law degree distribution, there are a few nodes with a very
high degree, which are too large to be stored in a block. For the exceptional
nodes, we divide and store them in multiple blocks.

In-memory layer. This layer is in charge of (1) loading blocks required for
the next operations into main memory, (2) aggregating the intermediate results
from the operation layer to generate the final result, and (3) writing the
final result to the storage layer. For this process, we define three buffers
and four index tables:

  * •

Input buffers ($B^{1}_{in}$ and $B^{2}_{in}$): storing rows of the input
matrices (i.e., blocks) required for the row-wise product.

  * •

Output buffer ($B_{out}$): storing the output result of a row-wise product
computed from the operation layer.

  * •

Buffer index tables ($T^{1}_{buf}$ and $T^{2}_{buf}$): indicating blocks
loaded in the two input buffers of $B^{1}_{in}$ and $B^{2}_{in}$.

  * •

Object index tables ($T^{1}_{obj}$ and $T^{2}_{obj}$): indicating the mappings
of blocks and objects in external storage.

The buffer index tables ($T^{*}_{buf}$) are referred to check whether the
required rows are already loaded on the main memory. If the rows are not on
the main memory, this layer requests ‘the block having the required rows’ to
the storage layer. In the object index tables ($T^{*}_{obj}$), the $i$-th
column represents the indices of objects stored in the $i$-th block. For
memory efficiency, we sort objects by their indices and store the only two
indices of the first and last objects. The object index tables ($T^{*}_{obj}$)
are pinned in the main memory to quickly load the necessary blocks. Note that
the pinned object index tables do not affect the entire performance of SAGE
since it requires only a small amount of memory less than $0.01\%$ of the
original graph.

Operation layer. This layer is in charge of performing actual multiplications
using the data in $B^{1}_{in}$ and $B^{2}_{in}$. This layer is performed as
follows: (1) reading two rows to be multiplied in $B^{1}_{in}$ and
$B^{2}_{in}$ by referring to $T^{1}_{buf}$ and $T^{2}_{buf}$, (2) performing a
row-wise product between them, and (3) writing the result (the intermediate
result of $M_{out}(i,;)$) to the output buffer. This process is repeated until
the input matrices are processed. The intermediate results are aggregated in
the in-memory layer to generate the final result $M_{\textit{out}}(i,;)$, and
then stored in the external storage.

Input : Input matrices $M^{1}_{in}$, $M^{2}_{in}$, Object index tables
$T^{1}_{obj}$, $T^{2}_{obj}$

Output : Output matrix $\mathcal{H}^{*}$

1 Function  _SAGE(_$M^{1}_{in}$ , $M^{2}_{in}$, $T^{1}_{obj}$,
$T^{2}_{obj}$_)_:

2 $T^{1}_{buf},T^{2}_{buf},M_{\textit{out}},\text{\sf
t\\_row}\leftarrow\emptyset$

3 for  _$M^{1}_{in}(i,;)\in M^{1}_{in}$_ do // processed in multiple threads

4 $\text{\sf LoadData}(i,T^{1}_{obj},T^{1}_{buf})$

5 for  _$M^{2}_{in}(j,;)\in M^{2}_{in}$_ do // $i$-th row computation

6 $\text{\sf LoadData}(j,T^{2}_{obj},T^{2}_{buf})$

7 $B_{out}\leftarrow B_{out}\cup\\{\text{\sf t\\_row}\\},\text{\sf
t\\_row}\leftarrow M^{\textit{1}}_{\textit{in}}(i,j)\times M^{2}_{in}(j,;)$

8

9 $M_{out}(i,;)\leftarrow\sum_{\text{\sf tmp}\in B_{out}}\text{\sf tmp}$ //
result aggregation

10 return _$M_{\textit{out}}$_

11Function  _LoadData(_$i,T_{obj},T_{buf}$_)_:

12 for  _$j=0,1,...,|T_{obj}|-1$_ do

13 $(n_{1},n_{2})\leftarrow T_{obj}[j]$ // first/last nodes in a block

14 if  _$n_{1}\leq i\leq n_{2}$_ then

15 $\text{\sf b\\_index}\leftarrow j$

16 break

17

18 if  _$\text{\sf b\\_index}\not\in T_{buf}$_ then // storage-memory I/O

19 $T_{buf}\leftarrow T_{buf}\cup\\{\text{\sf b\\_index}\\}$

20

21 return __

Algorithm 1 SAGE: Storage-based SpGEMM

####  3.2.2. Algorithm and performance consideration

Next, let us describe how SAGE performs SpGEMM based on the 3-layer
architecture. Algorithm 1 shows the entire process of SAGE (logical view).
Given the two input matrices, which are split and stored in multiple blocks,
and their object index tables, (Data loading) SAGE first loads the two rows to
be multiplied into the input buffers ($\text{\sf LoadData}(\cdot)$ in lines
10-18). Specifically, (1) the operation layer requests the two rows to the in-
memory layer, (2) the in-memory layer accesses the input buffer tables to
check whether the blocks having the required rows are already in the input
buffers, and (3) the in-memory layer sends the blocks having the required rows
to the operation layer if they are already in the input buffers, or otherwise,
requests the blocks to the storage layer. (Multiplication) Then, for the two
rows loaded in the main memory, SAGE computes the $i$-th row of the output
matrix, $M_{out}(i,;)$, by Eq. 1 (lines 3-8), where the intermediate results
(t_row) are stored in $B_{out}$ and aggregated to generate the final result
$M_{out}(i,;)$ in the in-memory layer. Finally, the in-memory layer sends
(i.e., writes) the final results of the output matrix to the storage layer.

SAGE employs a multi-thread based row-wise product to accelerate SpGEMM.
Basically, SAGE assigns each row of the former input matrix to each thread,
rather than the latter input matrix. Thus, the outer for loop in Algorithm 1
(lines 3-8) is performed by multiple threads in parallel. In this way, SAGE
can benefit from the parallelism of SpGEMM since the computation of each
thread is independent to each other. Furthermore, to maximize the parallelism,
we adopt a block-based workload allocation strategy to distribute workloads
into multiple threads evenly (elaborated in Section 3.3.1).

(a) Row-based allocation.

(b) Block-based allocation.

Figure 4. Comparison of (a) the row-based workload allocation and (b) the
block-based workload allocation.

###  3.3. Challenges and Strategies

We then point out three critical challenges leading to serious performance
degradation in SAGE and propose effective strategies to address them: (1)
block-based workload allocation, (2) in-memory partial aggregation, and (3)
distribution-aware memory allocation.

####  3.3.1. Workload balancing

To efficiently perform SpGEMM on large-scale graphs, existing works generally
adopt multi-thread based computation (Akbudak and Aykanat, 2017; Wang et al.,
2014; Xie et al., 2019; Nagasaka et al., 2018; Buluç and Gilbert, 2012; Azad
et al., 2016), where each thread is in charge of computing different rows
(nodes) of the output matrix in parallel (i.e., row-based workload
allocation). Thus, the workload of each thread is decided, depending on the
number of elements of each row. In real-world graphs, however, nodes (rows)
have a quite different number of edges (elements) as real-world graphs tend to
follow the power-law degree distribution. It implies that the workloads of
threads could be significantly different from each other. Such workload
differences among threads can adversely affect the overall performance of
SpGEMM by decreasing the benefit of parallel processing. Figure 4(a) shows an
example of row-base workload allocation for five rows with different lengths
(R1-R5) stored in five blocks (dotted boxes), where a large amount of
workloads for (R1), stored in three blocks, are assigned on a single thread
(T1) (i.e., over-loaded), while only a small amount of workloads for (R3, R4,
and R5), stored in a single block, are split and assigned on three threads
(T3-T5) (i.e., under-loaded). Thus, it is important to distribute workloads
evenly into threads to maximize the benefit of parallelizing SpGEMM (Challenge
1).

To tackle the challenge of workload balancing, we propose a simple yet
effective strategy that distributes workloads into multiple threads in a
block-based manner (block-based workload allocation). Specifically, SAGE
equipped with the block-based workload allocation assigns “the same number of
blocks” to each thread, rather than the same number of rows with different
numbers of edges, thereby distributing the workloads into threads evenly. Note
that this strategy is applied to the outer for loop (line 3) in Algorithm 1,
where the rows of the former input matrix stored in each block are assigned
together to the same thread. Figure 4(b) shows an example of our block-based
workload allocation, where the large row (R1), stored in three blocks, is
split and assigned on three threads (T1-T3), and the small rows (R3-R5),
grouped in a single block, are assigned together on a single thread (T5). As a
result, all five threads have the same amount of workloads (one block per
thread).

Theoretically, our block-based strategy improves the workload difference among
threads, compared to the row-based one, $O(|r|)\rightarrow O(b)$, where $|r|$
is the # of elements in a row $r$ and $b$ (fixed) is # of elements in a block
(in general, $|r|_{max}\gg b$).

Figure 5. The amount of storage-memory I/Os could be significantly reduced by
the in-memory partial aggregation.

####  3.3.2. Intermediate handling

As described in Section 3.2.2, the $i$-th row of the output matrix
($M_{out}(i,;)$) is computed by aggregating the intermediate results of
multiple row-wise products – the inner for loop (lines 5-8) in Algorithm 1.
For the aggregation, the intermediate results should be loaded on the output
buffer. The output buffer, however, may not be sufficient to store all
intermediates due to the limited size of the main memory, and it could be
overflowed by a large amount of intermediates. In that case, the intermediates
should be flushed into the storage multiple times to process the remaining
row-wise products. They are then loaded back into the main memory to be
aggregated in the end once the remaining row-wise products are completed. In
other words, a large amount of storage-memory I/Os could be generated by the
intermediates, which could degrade the performance of SpGEMM significantly.
Therefore, it is critical to handle a large amount of storage-memory I/Os
generated by the intermediates results (Challenge 2).

To address this challenge, we propose to column-wisely aggregate the
intermediate results of each row directly in the output buffer (in-memory
partial aggregation). Figure 5 illustrates the effect of our in-memory partial
aggregation strategy, where the intermediates with the same column indices are
aggregated in the main memory before they are flushed to storage. Accordingly,
SAGE is able to reduce a large amount of storage-memory I/Os unnecessarily
generated by the intermediates. Furthermore, we implement hash-based index
matching in SAGE to quickly find and aggregate the intermediates with the same
column index by following (Nagasaka et al., 2018). However, we have observed
that SpGEMM on extremely large-scale graphs may still cause even the partially
aggregated intermediates to be overflowed. In this case, SAGE stores them in
storage and aggregate them after finishing all of the row-wise products.

####  3.3.3. Memory management

SAGE manages three buffers in the main memory (i.e., two input and an output
buffers). Due to the limited size of the main memory, the amount of storage-
memory I/Os can vary depending on the proportions of the three buffers. For
example, (Case 1) if the input buffers are set too small, the data loaded into
the input buffers are processed in only a few iterations and replaced
frequently, incurring lots of storage-memory I/Os. On the other hand, (Case 2)
if the output buffer is set too small, it is easily overflowed by a small
number of intermediate results, which also causes storage-memory I/Os
frequently. Therefore, it is crucial to adjust the proportions of the three
buffers to reduce the amount of storage-memory I/Os (Challenge 3).

Toward this challenge, we propose a practical strategy that adjusts the
proportions of the three buffers, based on the characteristics of real-world
graphs (distribution-aware allocation). First, SAGE sets the size of the
former input buffer $B^{1}_{in}$ as ($b\times t$), where $b$ is the size of a
block and $t$ is the number of threads. Note that the size of ($b\times t$) is
sufficient for the former input buffer to load the blocks that all threads can
process simultaneously since SAGE allocates workloads to threads in a block-
based manner.

Second, SAGE determines the size of the output buffer $B_{out}$.
Unfortunately, it is infeasible to predict the exact size of the output row in
advance as it depends on data characteristics. Some previous works (Wang et
al., 2014; Xie et al., 2019; Akbudak and Aykanat, 2017; Nagasaka et al., 2018)
have tried to pre-compute the theoretical maximum size of the intermediate
results for each row, and set the output buffer size as the maximum size in
order to prevent the output buffer from being overflowed. However, it is
inappropriate for SpGEMM on real-world graphs following the power-law degree
distribution since the space required for the intermediate results of a
majority of rows is often much smaller than the theoretical maximum size
(i.e., inefficient memory usage). Thus, we argue that the output buffer with a
size much smaller than the maximum output size is practically sufficient to
process the intermediate results in most cases. Based on this intuition, we
can compute the maximum size of the intermediate results, $max\\_row$, and
sets the size of the output buffer as ($\alpha\times max\\_row\times t$),
considering the degree distribution of real-world graphs (i.e., distribution-
aware allocation), where $\alpha$ is the hyperparameter to adjust the output
buffer size.

Lastly, SAGE allocates the remaining portion of the main memory to the latter
input buffer $B^{2}_{in}$. Formally, given the capacity of the main memory
$C$, a block size $b$, the number of threads $t$, and the hyperparameter
$\alpha$, the main memory is allocated as:

(2) $\displaystyle{\textsf{SAGE}}\\_mem(C,b,t,\alpha)=\begin{cases}b\times
t,&B^{1}_{in}\\\ C-(|B^{1}_{in}|+|B_{out}|),&B^{2}_{in}\\\ \alpha\times
max\\_row\times t,&B_{out}.\end{cases}$

###  3.4. Complexity Analysis

The computational cost of SAGE derives mainly from (1) loading input matrices,
(2) row-wise products, (3) aggregating the intermediate results, and (4)
storing the final results in the output matrix. For brevity, we consider
SpGEMM between the two sparse matrices representing the same real-world graph
$G=(N,E)$, where $N$ and $E$ is a set of nodes and edges, respectively. The
input loading requires the time complexity of $O(|N|^{2})$ for two input
matrices. For $|N|$ many rows, the computational cost of each row-wise
products is $O(|N|\times\tilde{r})$, where $\tilde{r}$ is the average number
of elements in each row, where $\tilde{r}\ll|N|$ since a real-world graph
follows the power-law degree distribution. Thus, the overall complexity of
row-wise products is $O(|N|^{2}\times\tilde{r})$. The aggregation of the
intermediate results requires $O(|N|^{2}\times\tilde{r})$, i.e., (1) hash-
based index matching and (2) element-wise aggregation. Lastly, storing the
final results in the output matrix requires the complexity of
$O(|N|^{2}\times\tilde{r})$. As a result, the overall time complexity of SAGE
is $O(|N|^{2}\times\tilde{r})$, which is comparable to the time complexity of
common SpGEMM methods. We will empirically evaluate the scalability of SAGE
with the increasing size of graphs in Section 4.2.3.

##  4\. Experiments

In this section, we comprehensively evaluate SAGE by answering the following
evaluation questions (EQs):

  * •

EQ1: Does SAGE improve the performance of SpGEMM, compared to the existing
single-machine-based methods?

  * •

EQ2: Does SAGE improve the performance of SpGEMM, compared to the existing
distributed-system-based methods?

  * •

EQ3: How does the SpGEMM performance of SAGE scale up with the increasing size
of graphs?

  * •

EQ4: How effective are the proposed strategies of SAGE in improving the
performance of SpGEMM?

###  4.1. Experimental Setup

Table 3. Comparison of the single-machine-based approach with SAGE on 18 real-world datasets with various sizes. Dataset | Poisson3D | Enron | Epinions | Sphere | Filter3D | 598a | Torso2 | Cop20k | Cage12
---|---|---|---|---|---|---|---|---|---
Intel MKL | 0.1131 | 0.5792 | 0.6343 | 0.2521 | 0.3921 | 0.2660 | 0.0940 | 0.4120 | 0.3669
SAGE | 0.3206 | 0.7497 | 0.7001 | 0.5335 | 0.7436 | 0.5626 | 0.3976 | 0.7433 | 0.6331
Difference (sec.) |  0.2075$\uparrow$ |  0.1705$\uparrow$ |  0.0658$\uparrow$ |  0.2814$\uparrow$ |  0.3515$\uparrow$ |  0.2966$\uparrow$ |  0.3036$\uparrow$ |  0.3313$\uparrow$ |  0.2662$\uparrow$
Dataset | Slashdot | Gowalla | Pokec | Youtube | Livejournal | Wikipedia | UK-2005 | SK-2005 | Yahoo
Intel MKL | 1.9890 | 7.8299 | 12.8732 | N/A, OOM | N/A, OOM | N/A, OOM | N/A, OOM | N/A, OOM | N/A, OOM
SAGE | 1.7558 | 6.772 | 10.6391 | 38.734 | 94.8283 | 2,276 | 616 | 3,709 | 12,994
Difference (sec.) | 0.2332$\downarrow$ | 1.0579$\downarrow$ | 2.2341$\downarrow$ | - | - | - | - | - | -

Datasets. We use widely used 26 real-world datasets (Leskovec and Krevl, 2014;
Kolodziej et al., 2019) with various sizes, ranging from 10MB at the minimum
to 60GB at the maximum, and synthetic datasets. For the synthetic datasets, we
use the Recursive MATrix generator (R-MAT) (Chakrabarti et al., 2004) to
generate synthetic graphs following a power-law degree distribution. R-MAT
takes several parameters: (1) a benchmark type (Graph500, SSCA, or ER), (2)
the scale of the graph $n$, (3) the edge factor $e$, and (4) the skewness
parameters $a,b,c$, and $d$, determining the skewness of the generated graph.
Based on the input parameters, R-MAT generates a matrix
$M\in\mathbb{R}^{2^{n}\times 2^{n}}$ with $e\times 2^{n}$ elements. We use the
same parameters as used in (Demirci and Aykanat, 2020; Buluç and Gilbert,
2012; Azad et al., 2016).

Evaluation protocol. To comprehensively evaluate SAGE, we consider two
different types of SpGEMM: (1) SpGEMM between “the same two matrices”
($M\times M$) and (2) SpGEMM between “two different matrices” ($M_{1}\times
M_{2}$ or $M\times M^{T}$). The first type (1) is used in many graph
algorithms such as finding all-pairs shortest paths (D’alberto and Nicolau,
2007), self-similarity joins (He et al., 2010), and summarization of sparse
datasets (Ordonez et al., 2016), and the second type (2) is used in
collaborative filtering (Linden et al., 2003), triangle counting (Becchetti et
al., 2008), and similarity joins of two different sparse graphs (Ordonez et
al., 2016). For the evaluation metric, we use the execution time of SpGEMM of
each method because the goal of this work is to improve the performance of
graph algorithms by accelerating SpGEMM.

Competing methods. We compare SAGE with the following SpGEMM methods in our
experiments.

  * •

Intel MKL (Wang et al., 2014): the Intel MKL, a state-of-the-art single-
machine-based method, is the Intel CPU optimized open-source library
supporting a series of math functions such as SpGEMM and SpMV, which was also
used in (Xie et al., 2019; Akbudak and Aykanat, 2017; Nagasaka et al., 2018;
Srivastava et al., 2020; Pal et al., 2018).

  * •

SpSUMMA (Buluç and Gilbert, 2012): SpSUMMA is a distributed-system-based
method that parallelizes SpGEMM by partitioning input matrices into 2-D grids
and processing them using multiple machines.

  * •

Graphulo (Hutchison et al., 2015): Graphulo stores input matrices in a
distributed system exploiting the data locality to reduce the inter-machine
communication overhead.

  * •

gRRp (Demirci and Aykanat, 2020): gRRp adopts a bipartite graph-based
partitioning method for balancing workloads among multiple machines in a
distributed system.

Implementation details. We run all experiments on a single machine equipped
with Intel i7-7700K CPU, 64GB main memory, and 4TB SSD as the storage. We set
the number of threads $t$ as 4 to fully utilize all four physical cores in the
CPU. As explained in Section 3.3.3, SAGE allocates the main memory into three
buffers (i.e., two input and one output matrices) by the distributed-aware
allocation (Eq. 2), where we set the hyperparameter $\alpha$ as 12.5%. We will
empirically evaluate the impact of the hyperparameter $\alpha$ on the
performance of SAGE in Section 4.2.4.

###  4.2. Experimental Results

####  4.2.1. Single-machine-based approach (EQ1)

We first compare SAGE with the single-machine-based approach, Intel MKL (Wang
et al., 2014), in terms of the SpGEMM performance (i.e., execution time).

Setup. To comprehensively evaluate SAGE, we use 18 real-world datasets with
various sizes, 12 relatively small datasets that can be loaded on the main
memory and 6 large datasets that exceed the limit of the main memory capacity.
Then, we (1) perform SpGEMM using each method on the datasets for 5 times and
(2) measure the average execution time of SpGEMM of each method. Note that we
limit the main memory size to 16GB to rigorously evaluate the ability of SAGE
to handle the storage-memory I/Os.

Results and analysis. Table 3 shows the results. First, SAGE provides the
comparable performance to or even “better” than the Intel MKL, the in-memory
based method. Specifically, SAGE (1) finishes SpGEMM slightly slower (only
less than 0.4 sec.) than the Intel MKL in 9 datasets which are small enough to
load the entire data on the main memory and (2) even outperforms the Intel MKL
in Slashdot, Gowalla, and Pokec datasets. We highlight that these improvements
over the Intel MKL are significant, given that the Intel MKL is an in-memory
based SpGEMM method loading and processing the entire data on the main memory
(i.e., no storage-memory I/Os) and specially designed for the Intel CPU
architecture, while SAGE has to handle the inherent additional overhead for
data transfers between the main memory and external storage as described in
Section 3.2. Second, SAGE successfully performs SpGEMM on very large real-
world datasets (i.e., from Youtube to Yahoo datasets in Table 3) that the
Intel MKL could not handle due to the limited capacity of the main memory.

As a result, these results demonstrate that SAGE is able not only (1) to
perform SpGEMM on very large graphs by utilizing the sufficient capacity of
external storage (i.e., scalable), but also (2) to efficiently handle storage-
memory I/Os by addressing the three challenges, as described in Section 3.3.

Table 4. Statistics of synthetic datasets Benchmark |  Scale $n$ | # of nodes | # of edges
---|---|---|---
Graph500 (Graphulo and gRRp) | 15 | 32.7K | 524K
16 | 65.6K | 1.04M
17 | 131K | 2.09M
18 | 262K | 4.19M
SSCA (SpSUMMA) | 21 | 2.09M | 16.7M
22 | 4.19M | 33.5M
23 | 8.38M | 67.1M
24 | 16.7M | 134M
Table 5. Comparison of the distributed-system-based methods with SAGE on 13 real-world datasets with various sizes. Dataset | Sphere | Filter3D | 598a | Torso2 | Cage12 | 144 | Wave | Majorbasis | Scircuit | Mac | Offshore | Mario | Tmt_sym
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Graphulo | 7.37 | 12.91 | 9.57 | 4.81 | 13.28 | 11.97 | 11.85 | 8.07 | 6.36 | 4.63 | 18.86 | 9.07 | 13.55
gRRp | 1.50 | 3.71 | 1.65 | 0.86 | 2.63 | 2.35 | 1.87 | 1.52 | 1.02 | 1.31 | 3.89 | 1.70 | 3.34
SAGE | 0.36 | 0.77 | 0.34 | 0.11 | 0.54 | 0.48 | 0.41 | 0.25 | 0.18 | 0.25 | 0.93 | 0.30 | 0.62
Difference (sec.) | 1.14$\downarrow$ | 2.94$\downarrow$ | 1.31$\downarrow$ | 0.75$\downarrow$ | 2.09$\downarrow$ | 1.87$\downarrow$ | 1.46$\downarrow$ | 1.27$\downarrow$ | 0.84$\downarrow$ | 1.06$\downarrow$ | 2.96$\downarrow$ | 1.40$\downarrow$ | 2.72$\downarrow$

####  4.2.2. Distributed-system-based approach (EQ2)

In this experiment, we compare SAGE with three state-of-the-art distributed-
system-based SpGEMM methods: Graphulo (Hutchison et al., 2015), gRRp (Demirci
and Aykanat, 2020), and Sparse SUMMA (SpSUMMA) (Buluç and Gilbert, 2012). We
use their experimental results reported in (Demirci and Aykanat, 2020; Buluç
and Gilbert, 2012) by following (Kyrola et al., 2012; Roy et al., 2013; Jo et
al., 2019; Maass et al., 2017; Jun et al., 2018) since we use the exactly same
datasets as in (Demirci and Aykanat, 2020; Buluç and Gilbert, 2012) and it is
not feasible to construct the identical distributed systems that the existing
works used.

Setup. We use both real-world and synthetic datasets, which were used in
(Demirci and Aykanat, 2020; Buluç and Gilbert, 2012). Specifically, for
comparison with Graphulo and gRRp, we use 13 real-world and 4 synthetic
datasets, generated by Graph500 benchmark with the scale factor varying from
15 to 18, where we set the edge factor $e=16$ and the skewness parameters
$a=0.57$, $b=0.19$, $c=0.19$, and $d=0.05$ by following (Demirci and Aykanat,
2020). For comparison with SpSUMMA, we use 4 synthetic datasets, generated by
SSCA benchmark with the scale factor varying $n$ from 21 to 24, where we set
the edge factor $e=8$ and the skewness parameters $a=0.6$ and $b=c=d=0.4/3$ by
following (Buluç and Gilbert, 2012). Table 4 shows the statistics of the
synthetic datasets used in this experiment.

Then, we (1) perform SpGEMM using SAGE on the datasets for 5 times and (2)
measure the average execution time of SpGEMM of each method. We note that
Graphulo and gRRp have been conducted on the distributed system that consists
of 10 machines equipped with 2 Intel Xeon E5-2690 v4 processors connected via
the DGS-3120-24TC Ethernet (Demirci and Aykanat, 2020) and SpSUMMA has been
conducted on the Franklin XT4 system, where each machine is equipped with a
quad-core AMD Opteron processor and connected with each other via a 6.4 GB/s
network (Buluç and Gilbert, 2012).

Results and analysis. Table 5 shows the results of Graphulo, gRRp, and SAGE on
13 real-world datasets. The results reveal that SAGE consistently outperforms
the state-of-the-art distributed-system-based methods in all real-world
datasets. Specifically, SAGE finishes the SpGEMM by up to $7.7\times$ faster
and gRRp, the best performer among competitors. Whereas, Graphulo provides the
poor SpGEMM performance on real-world graphs (e.g., by up to $43\times$ slower
than SAGE in Torso2 dataset). This is because Graphulo, which adopts the outer
product, suffers from the significant overhead of indexing and aggregating a
large amount of the intermediate results of outer-product based SpGEMM as
described in Section 3.1. Although gRRp provides the performance of SpGEMM
much better than Graphulo, SAGE still significantly outperforms gRRp. We argue
that these performance improvement of SAGE over Graphulo and gRRp is
significant, considering they perform SpGEMM on the high-performance
distributed system that consists of 10 machines (Demirci and Aykanat, 2020).
Thus, these results imply that SAGE is able to efficiently handle large-scale
SpGEMM in a single machine, requiring only intra-machine communication
overhead (i.e., storage-memory I/Os), much smaller than the inter-machine
communication overhead.

GraphulogRRpSAGESpSUMMA-LSpSUMMA-MSpSUMMA-H

15161718$0$$5$$10$$15$$20$Graph scale factor ($n$)Execution time (sec.) (a)
Graph500

212223240204060Graph scale factor ($n$)Execution time (sec.) (b) SSCA.

Figure 6. Comparison of the distributed-system-based methods with SAGE on
synthetic datasets.

Figure 6 shows the results on synthetic datasets, generated by two benchmarks
(Graph500 and SSCA), where the $x$-axis represents the graph scale factor $n$
of each synthetic dataset and the $y$-axis represents the execution time
(sec.). Similar to the results on real-wold datasets, Figure 6(a) shows that
SAGE significantly outperforms the two distributed-system-based methods,
Graphulo and gRRp, across all graph scales. Moreover, as shown in Figure
6(b)111Buluç and Gilbert (2012) used three different number of cores in its
system with 3 (SpSUMMA-L), 4 (SpSUMMA-M), and 9 (SpSUMMA-H) machines,
respectively., SAGE successfully performs SpGEMM on very large-scale graphs
(the scale factor $n\geq 23$) that SpSUMMA-L and SpSUMMA-M fail to handle in a
single machine. As a result, These results demonstrate that SAGE is able to
process large-scale real-world graphs successfully with only limited computing
resources of a single machine, rather than those of costly distributed-systems
(i.e., cost-effective).

HighNormalLow

182022242628$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$Graph scale
factorExecution time (sec.) (a) Graph500

182022242628$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$Graph scale
factorExecution time (sec.) (b) ER

Figure 7. The SpGEMM execution time of SAGE with the increasing sizes of
synthetic graphs.

####  4.2.3. Scalability (EQ3).

In this experiment, we evaluate the scalability of SAGE with the increasing
sizes of graphs.

Setup. We use synthetic datasets, generated by two different benchmarks
(Graph500 and ER) (Azad et al., 2016) with various graph scales $n$ ranging
from 18 to 28, where the scale factor $n$ decides the number of nodes in a
generated graph (i.e., $|N|=2^{n}$). For each generated graph, we also
consider two different edge factors to control the number of edges (i.e., one
for a high-degree graph and the other for a low-degree graph). Specifically,
we set the edge factor for the high-degree and low-degree graphs as $1.5\times
e$ and $0.5\times e$, respectively222We could not show the results of high-
degree graphs with graph scale 28 since R-MAT generator could not generate
them due to the integer overflow.. Then, we perform SpGEMM (type (2)) on the
generated graphs for 5 times and measure the average execution time (sec.).

Results and analysis. Figure 7 shows the results, where the $x$-axis
represents the scale factor $n$, the $y$-axis represents the execution time
(sec. in log-scale), and the red dotted vertical line represents the limits of
existing single-machine-based methods. We observe that the execution time of
SAGE tends to increase by 4.5 times as the scale factor increases by 2 (i.e.,
a graph gets $4\times$ larger). This result demonstrates that the SpGEMM
performance of SAGE successfully scales up to very large-scale graphs,
exceeding the capacity of the main memory (red dotted vertical line), and SAGE
provides (almost) linear scalability with the increasing sizes of graphs.

WikipediaUK-2005SK-2005Yahoo05k10k15kExecution time(sec.)SAGE-NoSAGE-BWSAGE-
PASAGE-AllO.O.S.O.O.S.O.O.S.O.O.S. (a) Effects on the SpGEMM performance.

WikipediaUK-2005SK-2005Yahoo050100Relative amountof I/Os (%)SAGE w/o PASAGE w/
PA (b) Analysis of the in-memory-partial aggregation.

Figure 8. Effects of our proposed strategies on SAGE.

####  4.2.4. Ablation study (EQ4).

Lastly, in this section, we evaluate the effectiveness of the proposed
strategies of SAGE individually.

Setup (1). We first evaluate the first two strategies: (1) block-based
workload allocation and (2) in-memory partial aggregation. To rigorously
evaluate the two strategies, we (1) use four largest real-world datasets
(i.e., Wikipedia, UK-2005, SK-2005, and Yahoo datasets (Jo et al., 2019)),
which could not be loaded into main memory, and (2) limit the size of the main
memory to 2GB to incur a large amount of storage-memory I/Os. We compare the
following four versions of SAGE:

  * •

SAGE-No: a baseline without any optimizations

  * •

SAGE-BW: the version with the block-based workload allocation

  * •

SAGE-PA: the version with the in-memory partial aggregation

  * •

SAGE-All: the final version with both optimizations

We perform SpGEMM (type (1)) using each version and measure the execution time
(sec.). For this experiment, we set the hyperparameter $\alpha=12.5\%$ for the
main memory buffer allocation.

Results and analysis (1). As shown in Figure 8(a), each of the proposed
strategy is effective in improving the performance of SpGEMM and SAGE-All
consistently provides the best performance in all datasets (up to $4.7\times$
over SAGE-no). Also, for the two largest datasets (i.e., SK-2005 and Yahoo),
SAGE-no and SAGE-BW fail to perform SpGEMM because the size of intermediate
results even exceeds the capacity of the 4TB external storage, where ‘O.O.S.’
indicates ‘out of storage’. Interestingly, SAGE-PA consistently outperforms
SAGE-BW in all cases, which implies that, in case of large-scale SpGEMM,
handling storage-memory I/Os could be more critical than balancing workloads,
in particular when the main memory size of a single machine is limited.

For more analysis of the in-memory partial aggregation, we measure the amount
of storage-memory I/Os generated by the intermediate results during SpGEMM of
SAGE with/without the in-memory partial aggregation (SAGE w/ PA and SAGE w/o
PA in Figure 8(b)). Figure 8(b) shows the relative amount of storage-memory
I/Os, where our in-memory partial aggregation significantly reduces the amount
of storage-memory I/Os (up to $87\%$ in UK-2005). This result demonstrates
that in performing storage-based SpGEMM on large-scale graphs, a substantial
amount of storage-memory I/Os could be generated; thus, it is crucial to
efficiently handle these I/Os in order to achieve good performance, as we
claimed in Section 3.3.

Setup (2). We also evaluate our distribution-aware memory allocation. As
described in in Section 3.3.3, we set (1) the former input buffer
$B^{1}_{in}=b\times t$, where $b$ is the size of a block (1MB) and $t$ is the
number of threads; and (2) the output buffer $B_{out}=\alpha\times
max\\_row\times t$, where $max\\_row$ here indicates the maximum size of the
final results for each row. Then, we perform SpGEMM with varying (1) the
numbers of threads (2, 4, and 8), (2) the size of $B^{1}_{in}$ (1MB to 16MB),
and (3) the hyperparameter $\alpha$ (6.25% to 50%).

Results and analysis (2). Figure 9(a) shows that the performance of SAGE gets
improved as the number of threads and the size of $B^{1}_{in}$ increase, where
the performance improvement converges at a specific point where the size of
$B^{1}_{in}$ is equal to $b\times t$, which implies that the buffer size is
sufficient to load the blocks for that all threads can process simultaneously,
as explained in Section 3.3.3. Figure 9(b) shows the performance of SAGE with
respect to the size of main memory and the hyperparameter $\alpha$. The
performance of SAGE tends to be improved as $\alpha$ and the size of the main
memory increase, but the performance improvement gets smaller as $\alpha$
increases, which implies that about $12.5\sim 25\%$ of the maximum size is
enough for $B_{out}$ to process the intermediate results of a majority of
rows. As a result, we believe that our distribution-aware memory allocation is
effective and practical in large-scale SpGEMM.

$1$$2$$4$$8$$16$$0$$200$$400$the size of $B^{1}_{in}$ (MB)Execution time
(sec.)t=2t=4t=8 (a) Input buffer

$6.25$$12.5$$25$$50$2.2k2.4k2.6khyperparameter $\alpha$ (%)0.5GB1GB2GB4GB (b)
Output buffer

Figure 9. Effects of distribution-aware memory allocation.

##  5\. Conclusions

In this paper, we point out the limitations of previous approaches for SpGEMM:
(i) the single-machine-based approach cannot handle large-scale graphs
surpassing the capacity of the main memory (i.e., not scalable), and (ii) the
distributed-system-based approach requires a substantial amount of inter-
machine communication overhead (i.e., not efficient). To address both
challenges of scalability and efficiency, we propose a novel storage-based
approach to SpGEMM (SAGE) with the 3-layer architecture to efficiently handle
storage-memory I/Os. We further identify three important challenges that could
cause serious performance degradation in storage-based SpGEMM and propose
three effective strategies to address them: (1) block-based workload
allocation, (2) in-memory partial aggregation, and (3) distribution-aware
memory allocation. Through comprehensive evaluation, we demonstrate the
superiority of SAGE in terms of (1) scalability, (2) efficiency, and (3)
effectiveness.

##  6\. Acknowledgments

This is a joint work between Samsung Electronics Co., Ltd and Hanyang
University. This work was supported by Institute of Information &
communications Technology Planning & Evaluation (IITP) grant funded by the
Korea government (MSIT) (No.2022-0-00352 and No.RS-2022-00155586).

## References

  * (1)
  * Akbudak and Aykanat (2017) Kadir Akbudak and Cevdet Aykanat. 2017.  Exploiting locality in sparse matrix-matrix multiplication on many-core architectures.  _IEEE Transactions on Parallel and Distributed Systems_ 28, 8 (2017), 2258–2271.
  * Anderson et al. (2016) Michael J Anderson, Narayanan Sundaram, Nadathur Satish, Md Mostofa Ali Patwary, Theodore L Willke, and Pradeep Dubey. 2016\.  Graphpad: Optimized graph primitives for parallel and distributed platforms. In _2016 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 313–322.
  * Azad et al. (2016) Ariful Azad, Grey Ballard, Aydin Buluc, James Demmel, Laura Grigori, Oded Schwartz, Sivan Toledo, and Samuel Williams. 2016\.  Exploiting multiple levels of parallelism in sparse matrix-matrix multiplication.  _SIAM Journal on Scientific Computing_ 38, 6 (2016), C624–C651.
  * Baek et al. (2021) Daehyeon Baek, Soojin Hwang, Taekyung Heo, Daehoon Kim, and Jaehyuk Huh. 2021.  InnerSP: A Memory Efficient Sparse Matrix Multiplication Accelerator with Locality-Aware Inner Product Processing. In _2021 30th International Conference on Parallel Architectures and Compilation Techniques (PACT)_. IEEE, 116–128.
  * Becchetti et al. (2008) Luca Becchetti, Paolo Boldi, Carlos Castillo, and Aristides Gionis. 2008.  Efficient semi-streaming algorithms for local triangle counting in massive graphs. In _Proceedings of the ACM International Conference on Knowledge Discovery and Data Mining_. ACM, 16–24.
  * Buluç and Gilbert (2012) Aydin Buluç and John R Gilbert. 2012.  Parallel sparse matrix-matrix multiplication and indexing: Implementation and experiments.  _SIAM Journal on Scientific Computing_ 34, 4 (2012), C170–C191.
  * Chakrabarti et al. (2004) Deepayan Chakrabarti, Yiping Zhan, and Christos Faloutsos. 2004\.  R-MAT: A recursive model for graph mining. In _Proceedings of the 2004 SIAM International Conference on Data Mining_. SIAM, 442–446.
  * D’alberto and Nicolau (2007) Paolo D’alberto and Alexandru Nicolau. 2007.  R-Kleene: A high-performance divide-and-conquer algorithm for the all-pair shortest path for densely connected networks.  _Algorithmica_ 47, 2 (2007), 203–213.
  * Davis (2019) Timothy A Davis. 2019\.  Algorithm 1000: SuiteSparse: GraphBLAS: Graph algorithms in the language of sparse linear algebra.  _ACM Transactions on Mathematical Software (TOMS)_ 45, 4 (2019), 1–25.
  * Demirci and Aykanat (2020) Gunduz Vehbi Demirci and Cevdet Aykanat. 2020.  Scaling sparse matrix-matrix multiplication in the accumulo database.  _Distributed and Parallel Databases_ 38, 1 (2020), 31–62.
  * Deveci et al. (2017) Mehmet Deveci, Christian Trott, and Sivasankaran Rajamanickam. 2017\.  Performance-portable sparse matrix-matrix multiplication for many-core architectures. In _2017 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)_. IEEE, 693–702.
  * Gremse et al. (2018) Felix Gremse, Kerstin Kupper, and Uwe Naumann. 2018\.  Memory-efficient sparse matrix-matrix multiplication by row merging on many-core architectures.  _SIAM Journal on Scientific Computing_ 40, 4 (2018), C429–C449.
  * Gu et al. (2017) Rong Gu, Yun Tang, Chen Tian, Hucheng Zhou, Guanru Li, Xudong Zheng, and Yihua Huang. 2017.  Improving execution concurrency of large-scale matrix multiplication on distributed data-parallel platforms.  _IEEE Transactions on Parallel and Distributed Systems_ 28, 9 (2017), 2539–2552.
  * Hamedani and Kim (2017) Masoud Reyhani Hamedani and Sang-Wook Kim. 2017.  JacSim: An accurate and efficient link-based similarity measure in graphs.  _Information Sciences_ 414 (2017), 203–224.
  * Han et al. (2013) Wook-Shin Han, Sangyeon Lee, Kyungyeol Park, Jeong-Hoon Lee, Min-Soo Kim, Jinha Kim, and Hwanjo Yu. 2013.  TurboGraph: A fast parallel graph engine handling billion-scale graphs in a single PC. In _Proceedings of the ACM international conference on knowledge discovery and data mining (KDD)_. 77–85.
  * He et al. (2010) Guoming He, Haijun Feng, Cuiping Li, and Hong Chen. 2010\.  Parallel SimRank computation on large graphs with iterative aggregation. In _Proceedings of the ACM international conference on knowledge discovery and data mining (SIGKDD)_. 543–552.
  * Hutchison et al. (2015) Dylan Hutchison, Jeremy Kepner, Vijay Gadepally, and Adam Fuchs. 2015.  Graphulo implementation of server-side sparse matrix multiply in the Accumulo database. In _2015 IEEE High Performance Extreme Computing Conference (HPEC)_. IEEE, 1–7.
  * Jo et al. (2019) Yong-Yeon Jo, Myung-Hwan Jang, Sang-Wook Kim, and Sunju Park. 2019\.  RealGraph: A graph engine leveraging the power-law distribution of real-world graphs. In _Proceedings of the World Wide Web Conference_. ACM, 807–817.
  * Jun et al. (2018) Sang-Woo Jun, Andy Wright, Sizhuo Zhang, Shuotao Xu, et al. 2018\.  GraFBoost: Using accelerated flash storage for external graph analytics. In _2018 ACM/IEEE 45th Annual International Symposium on Computer Architecture (ISCA)_. IEEE, 411–424.
  * Kang et al. (2020) Yoonsuk Kang, Jun Seok Lee, Won-Yong Shin, and Sang-Wook Kim. 2020.  Cr-graph: Community reinforcement for accurate community detection. In _Proceedings of the 29th ACM International Conference on Information & Knowledge Management_. 2077–2080.
  * Karypis and Kumar (1998) George Karypis and Vipin Kumar. 1998.  Multilevel k-way partitioning scheme for irregular graphs.  _J. Parallel and Distrib. Comput._ 48, 1 (1998), 96–129.
  * Kepner et al. (2016) Jeremy Kepner, Peter Aaltonen, David Bader, Aydin Buluç, Franz Franchetti, John Gilbert, Dylan Hutchison, Manoj Kumar, Andrew Lumsdaine, Henning Meyerhenke, et al. 2016\.  Mathematical foundations of the GraphBLAS. In _2016 IEEE High Performance Extreme Computing Conference (HPEC)_. IEEE, 1–9.
  * Kepner and Gilbert (2011) Jeremy Kepner and John Gilbert. 2011.  _Graph algorithms in the language of linear algebra_.  SIAM.
  * Kleinberg (1999) Jon M Kleinberg. 1999\.  Authoritative sources in a hyperlinked environment.  _J. ACM_ 46, 5 (1999), 604–632.
  * Ko et al. (2021a) Yunyong Ko, Kibong Choi, Jiwon Seo, and Sang-Wook Kim. 2021a.  An In-Depth Analysis of Distributed Training of Deep Neural Networks. In _Proceedings of the IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 994–1003.
  * Ko et al. (2021b) Yunyong Ko, Jae-Seo Yu, Hong-Kyun Bae, Yongjun Park, Dongwon Lee, and Sang-Wook Kim. 2021b.  MASCOT: A Quantization Framework for Efficient Matrix Factorization in Recommender Systems. In _2021 IEEE International Conference on Data Mining (ICDM)_. IEEE, 290–299.
  * Kolodziej et al. (2019) Scott P Kolodziej, Mohsen Aznaveh, Matthew Bullock, Jarrett David, Timothy A Davis, Matthew Henderson, Yifan Hu, and Read Sandstrom. 2019\.  The suitesparse matrix collection website interface.  _Journal of Open Source Software_ 4, 35 (2019), 1244\.
  * Kyrola et al. (2012) Aapo Kyrola, Guy E Blelloch, and Carlos Guestrin. 2012\.  GraphChi: Large-scale graph computation on just a pc. In _Proceedings of the USENIX symposium on operating systems design and implementation (OSDI)_. 31–46.
  * Leskovec et al. (2010a) Jure Leskovec, Daniel Huttenlocher, and Jon Kleinberg. 2010a.  Predicting positive and negative links in online social networks. In _Proceedings of the 19th international conference on World wide web_. 641–650.
  * Leskovec et al. (2010b) Jure Leskovec, Daniel Huttenlocher, and Jon Kleinberg. 2010b.  Signed networks in social media. In _Proceedings of the SIGCHI conference on human factors in computing systems_. 1361–1370.
  * Leskovec et al. (2007) Jure Leskovec, Jon Kleinberg, and Christos Faloutsos. 2007\.  Graph evolution: Densification and shrinking diameters.  _ACM Transactions on Knowledge Discovery from Data_ 1, 1 (2007), 1–41.
  * Leskovec and Krevl (2014) Jure Leskovec and Andrej Krevl. 2014.  SNAP Datasets: Stanford Large Network Dataset Collection.  http://snap.stanford.edu/data.
  * Linden et al. (2003) Greg Linden, Brent Smith, and Jeremy York. 2003.  Amazon. com recommendations: Item-to-item collaborative filtering.  _IEEE Internet computing_ 7, 1 (2003), 76–80.
  * Liu et al. (2016) Hang Liu, H Howie Huang, and Yang Hu. 2016.  ibfs: Concurrent breadth-first search on gpus. In _Proceedings of the 2016 International Conference on Management of Data_. 403–416.
  * Maass et al. (2017) Steffen Maass, Changwoo Min, Sanidhya Kashyap, Woonhak Kang, Mohan Kumar, and Taesoo Kim. 2017\.  Mosaic: Processing a trillion-edge graph on a single machine. In _Proceedings of the Twelfth European Conference on Computer Systems_. 527–543.
  * Nagasaka et al. (2018) Yusuke Nagasaka, Satoshi Matsuoka, Ariful Azad, and Aydin Buluç. 2018.  High-performance sparse matrix-matrix products on intel knl and multicore architectures. In _Proceedings of the 47th International Conference on Parallel Processing Companion_. 1–10.
  * Nagasaka et al. (2017) Yusuke Nagasaka, Akira Nukada, and Satoshi Matsuoka. 2017\.  High-performance and memory-saving sparse general matrix-matrix multiplication for nvidia pascal gpu. In _2017 46th International Conference on Parallel Processing (ICPP)_. IEEE, 101–110.
  * Ordonez et al. (2016) Carlos Ordonez, Yiqun Zhang, and Wellington Cabrera. 2016\.  The Gamma matrix to summarize dense and sparse data sets for big data analytics.  _IEEE Transactions on Knowledge and Data Engineering_ 28, 7 (2016), 1905–1918.
  * Page et al. (1999) Lawrence Page, Sergey Brin, Rajeev Motwani, and Terry Winograd. 1999.  _The PageRank citation ranking: Bringing order to the web_.  Technical Report. Stanford InfoLab.
  * Pal et al. (2018) Subhankar Pal, Jonathan Beaumont, Dong-Hyeon Park, Aporva Amarnath, Siying Feng, Chaitali Chakrabarti, Hun-Seok Kim, David Blaauw, Trevor Mudge, and Ronald Dreslinski. 2018\.  Outerspace: An outer product based sparse matrix multiplication accelerator. In _2018 IEEE International Symposium on High Performance Computer Architecture (HPCA)_. IEEE, 724–736.
  * Roy et al. (2013) Amitabha Roy, Ivo Mihailovic, and Willy Zwaenepoel. 2013\.  X-Stream: Edge-centric graph processing using streaming partitions. In _Proceedings of the ACM symposium on operating systems principles (SOSP)_. 472–488.
  * Sedgewick and Wayne (2011) Robert Sedgewick and Kevin Wayne. 2011.  _Algorithms_.  Addison-Wesley professional.
  * Selvitopi et al. (2019) Oguz Selvitopi, Gunduz Vehbi Demirci, Ata Turk, and Cevdet Aykanat. 2019.  Locality-aware and load-balanced static task scheduling for MapReduce.  _Future Generation Computer Systems_ 90 (2019), 49–61.
  * Srivastava et al. (2020) Nitish Srivastava, Hanchen Jin, Jie Liu, David Albonesi, and Zhiru Zhang. 2020.  Matraptor: A sparse-sparse matrix multiplication accelerator based on row-wise product. In _2020 53rd Annual IEEE/ACM International Symposium on Microarchitecture (MICRO)_. IEEE, 766–780.
  * Suzuki et al. (2003) Kenji Suzuki, Isao Horiba, and Noboru Sugie. 2003\.  Linear-time connected-component labeling based on sequential local operations.  _Computer Vision and Image Understanding_ 89, 1 (2003), 1–23.
  * Wang et al. (2014) Endong Wang, Qing Zhang, Bo Shen, Guangyong Zhang, Xiaowei Lu, Qing Wu, and Yajuan Wang. 2014.  Intel math kernel library.  In _High-Performance Computing on the Intel® Xeon Phi™_. Springer, 167–188.
  * Xie et al. (2019) Zhen Xie, Guangming Tan, Weifeng Liu, and Ninghui Sun. 2019\.  IA-SpGEMM: An input-aware auto-tuning framework for parallel sparse matrix-matrix multiplication. In _Proceedings of the ACM International Conference on Supercomputing_. 94–105.
  * Zhang et al. (2020) Zhekai Zhang, Hanrui Wang, Song Han, and William J Dally. 2020\.  Sparch: Efficient architecture for sparse matrix multiplication. In _2020 IEEE International Symposium on High Performance Computer Architecture (HPCA)_. IEEE, 261–274.
  * Zheng et al. (2015) Da Zheng, Disa Mhembere, Randal Burns, Joshua Vogelstein, Carey E Priebe, and Alexander S Szalay. 2015.  FlashGraph: Processing billion-node graphs on an array of commodity SSDs. In _Proceedings of the USENIX conference on file and storage technologies (FAST)_. 45–58.
  * Zheng et al. (2016) Da Zheng, Disa Mhembere, Vince Lyzinski, Joshua T Vogelstein, Carey E Priebe, and Randal Burns. 2016\.  Semi-external memory sparse matrix multiplication for billion-node graphs.  _IEEE Transactions on Parallel and Distributed Systems_ 28, 5 (2016), 1470–1483.


# Exact Renormalization of a Scalar Field

Kang-Sin Choi<EMAIL_ADDRESS>Scranton Honors Program, Ewha Womans
University, Seoul 03760, Korea  Institute of Mathematical Sciences, Ewha
Womans University, Seoul 03760, Korea

###### Abstract

Using Wilsonian renormalization, we calculate the quantum correction to
observable quantities, not bare parameters, of a scalar field. A physical
parameter (a mass-squared or a quartic coupling) at an energy scale $\mu$ is
obtained from that at a reference scale by integrating in the degrees of
freedom in between. Heavy modes decouple and the ultraviolet scale dependence
is canceled in the observables. An exact renormalization group equation is
parametrized by the low-energy scale $\mu$.

Renormalization in quantum field theory informs us how to deal with the energy
dependence of physical parameters. In computing their quantum corrections,
divergences arise. They have become under control by Wilsonian renormalization
Wilson:1971bg ; Polchinski:1983gv . The related cutoff is interpreted as the
energy scale or the dimensionful parameter of new physics. We can follow how
such parameters inherit to the low-energy ones.

So far, the renormalization procedure has focused on the bare parameters of a
given quantum field theory Choi:2023cqs . In this letter, we perform Wilsonian
renormalization for physically observable parameters, especially the scalar
masses. For instance, if we take a reference as the pole mass, its quantum
correction is expressed as the difference of the bare parameters at different
energy scales, canceling the cutoff dependence Choi:2023cqs .

Take a Higgs scalar $\phi(x)$ couped to a heavy scalar $X(x)$ and a fermion
$t(x)$. We consider a Euclidian action

$\begin{split}-{\cal
L}&\supset\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi+\frac{1}{2}m_{0}^{2}\phi^{2}+\frac{\lambda_{0}}{4}\phi^{4}\\\
&\quad+\frac{1}{2}\partial_{\mu}X\partial^{\mu}X+\frac{1}{2}M_{0}^{2}X^{2}+\frac{\kappa_{0}}{4}\phi^{2}X^{2}\\\
&\quad+\overline{t}_{L,R}\gamma^{\mu}\partial_{\mu}t_{L,R}+\frac{y_{t,0}}{\sqrt{2}}(v_{0}+\phi)(\overline{t}_{L}t_{R}+\overline{t}_{R}t_{L}).\end{split}$
(1)

Here, we expanded the Higgs field around the electroweak vacuum expectation
value (VEV) $v_{0}$. Thus it is a real scalar field usually denoted by $h$.
Other terms are omitted for simplicity. We consider a hierarchy of energy
scales

$\Lambda<M_{0}<\Lambda_{?},$ (2)

that is, we assume that the above Lagrangian is valid below an energy scale
$\Lambda_{?}$. We shall be interested in the scale way below the physics of
the field $X$. In principle, we do not need $\Lambda$, but in practice, we do
not know the necessary dimensionful parameter $M_{0}$. For $\phi$, we consider
both cases $m_{0}\ll M_{0}$ and $m_{0}\sim M_{0}$.111We should deal with the
renormalized masses, not the bare ones: see below. We are interested in the
mass and the quartic coupling of the scalar field $\phi$, so except for those,
we assume that the renormalized couplings are obtained in the same way. The
corresponding quantities are denoted with the subscript as the scale.

We calculate the mass of the field $\phi$ in the low energy. Since the mass
$M$ of the field $X$ is very large from the low-energy viewpoint, we expect it
to be decoupled. We also separate, in the scalar field, the high-frequency
part $\hat{\phi}(k)$, whose Euclidianized momentum is larger than $\Lambda$,
as

$\phi(k)\to\phi(k)+\hat{\phi}(k),$ (3)

and let the remaining one $\phi(k)$ fluctuate below $\Lambda$ Wilson:1971bg .
The relevant interactions are

$V=\frac{3\lambda_{0}}{2}\phi^{2}\hat{\phi}^{2}+\frac{\lambda_{0}}{4}\hat{\phi}^{4}+\frac{\kappa_{0}}{4}\phi^{2}X^{2}+\frac{y_{t,0}}{\sqrt{2}}\phi(\overline{t}_{L}t_{R}+\overline{t}_{R}t_{L}).$
(4)

The effective mass operator is obtained by contracting two $\hat{\phi}$ fields
in the interactions (4), using the Feynman propagator Polchinski:1983gv . We
do not consider external momenta. For the high-frequency scalar, we have

$\begin{split}\Sigma_{\hat{\phi}}(\mu^{2})&=\frac{3}{2}\int_{\mu<|k|<\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\frac{\lambda_{k}}{k^{2}+m_{k}^{2}}\\\
&=\frac{3}{2}\int_{\mu^{2}}^{\Lambda^{2}}\frac{dk^{2}k^{2}}{(4\pi)^{2}}\frac{\lambda_{k}}{k^{2}+m_{k}^{2}}.\\\
\end{split}$ (5)

Here, we took an arbitrary low energy scale $\mu<\Lambda$ and used the
nontrivial propagator only in the range indicated in the integral. In
Wilsonian renormalization, the cutoff $\Lambda$ is not the regularized
infinity but the scale we specify. In the propagator, we used a renormalized
mass $m_{k}$ at the energy scale $k$ and a quartic coupling $\lambda_{k}$ that
we clarify shortly. This makes the definition of mass recursive, but we can
approximate and calculate it perturbatively later.

Similarly, the heavy field $X$ and the top quark contribute to the mass by
integrating out over the same region

$\displaystyle\tilde{\Sigma}_{X}(\mu^{2})$
$\displaystyle=\frac{1}{4}\int_{\mu^{2}}^{\Lambda^{2}}\frac{dk^{2}k^{2}}{(4\pi)^{2}}\frac{\kappa_{k}}{k^{2}+M_{k}^{2}},$
(6) $\displaystyle\tilde{\Sigma}_{t}(\mu^{2})$
$\displaystyle=-\int_{\mu^{2}}^{\Lambda^{2}}\frac{dk^{2}k^{2}}{(4\pi)^{2}}\frac{y_{t}^{2}(k^{2}-m_{t,k}^{2})}{(k^{2}+m_{t,k}^{2})^{2}}.$
(7)

In the low-energy theory, the mass is corrected as

$m_{\mu}^{2}\equiv m_{0}^{2}+\sum\tilde{\Sigma}_{i}(\mu^{2}),$ (8)

where the summation runs over all the contributions we consider.

Apparently, the mass correction depends on the high scale $M$ and the cutoff
$\Lambda$ quadratically. We can further infer that all the bare parameters,
including $m_{0}^{2},\lambda_{0}$ must also be dependent on the cutoff
$\Lambda_{?}$ and unknown physics beyond that. Therefore, the mass is not
well-defined in low energy below $\Lambda$. This is the hierarchy problem.

However, this has been about the bare parameters. In what follows, we show
that we can mention observable mass and compute the quantum correction of it,
whose result does not depend on high energy Choi:2023cqs .

Now we extract what we can observe. Only the combination (8) can be
observable, so firstly we express it in reference to the pole mass. The
Feynman propagator for the low-energy scalar $\phi(x)$ with the momentum $k$
has a pole at the mass $k^{2}=m_{k}^{2}$ defined in (8). We define a “pole
mass” $m_{h}$ as that satisfying

$m_{h}^{2}=m_{0}^{2}+\sum\tilde{\Sigma}_{i}(m_{h}^{2}).$ (9)

This is a natural reference point, but any reference would be good. Since here
we do not consider the kinematics of the field, the mass (8) is independent of
$k$ and the field renormalization is not necessary.

We can express the effective mass at scale $\mu$ in terms of the pole mass
Choi:2023cqs

$\begin{split}m_{\mu}^{2}&=m_{h}^{2}+\sum\left[\tilde{\Sigma}_{i}(\mu^{2})-\tilde{\Sigma}_{i}(m_{h}^{2})\right]\\\
&\equiv m_{h}^{2}+\delta m_{h}^{2}(\mu^{2}).\end{split}$ (10)

The unobservable bare mass $m_{0}$ is removed and the mass is understood as
the quantum correction from the pole mass. Note that it depends on the scale
$\mu$ we consider: like the upper limit $\Lambda$, we specified to which
energy $\mu$ we run down. We consider quantum correction from this reference
mass we can measure by experiment. This reference plays a similar role as the
renormalization condition. The mass at different energy scales is now physical
in the sense that we can measure it as well.

Now, the mass correction consists of the combinations. For the high-frequency
scalar $\hat{\phi}$, it is

$\tilde{\Sigma}_{\hat{\phi}}(\mu^{2})-\tilde{\Sigma}_{\hat{\phi}}(m_{h}^{2})=\frac{3}{2}\int^{m_{h}^{2}}_{\mu^{2}}\frac{dk^{2}k^{2}}{(4\pi)^{2}}\frac{\lambda_{k}}{k^{2}+m_{k}^{2}}.$
(11)

This has a natural interpretation, faithfully following the Wilsonian program:
The physical parameter (the mass squared) at scale $\mu^{2}$ is obtained from
that at $m_{h}^{2}$ by integrating in the degrees of freedom from $m_{h}^{2}$
up to $\mu^{2}$. Considering multiple parameters, the renormalization group
flow may branch, so we should technically understand this as the minus of the
integrating out. Obviously, the integration is finite, as if we did not need
any regularization.

Since the running interval in the energy scale is short, we can approximate
the mass to be constant around the Higgs pole mass. Also, we show shortly that
the quartic coupling does not run much, so we approximate the dimensionless
parameters as constants. Then, the mass correction (10) becomes

$\begin{split}\delta
m(\mu^{2})&\simeq\frac{\kappa}{64\pi^{2}}\left[-\mu^{2}+m_{h}^{2}-M^{2}\log\frac{m_{h}^{2}+M^{2}}{\mu^{2}+M^{2}}\right]\\\
&+\frac{3\lambda}{32\pi^{2}}\left[-\mu^{2}+m_{h}^{2}-m_{h}^{2}\log\frac{2m_{h}^{2}}{\mu^{2}+m_{h}^{2}}\right]\\\
&+\frac{y_{t}^{2}}{16\pi^{2}}\bigg{[}-\mu^{2}+m_{h}^{2}-3m_{t}^{2}\log\frac{\mu^{2}+m_{t}^{2}}{m_{h}^{2}+m_{t}^{2}}\\\
&\quad-2m_{t}^{4}\left.\left(\frac{1}{m_{h}^{2}+m_{t}^{2}}-\frac{1}{\mu^{2}+m_{t}^{2}}\right)\right],\end{split}$
(12)

All the couplings are matched at $m_{h}$ and written down without subscripts.
The dependences on $\Lambda$, the quadratic and the logarithmic, are removed,
as promised.

Figure 1: Decoupling in the correction of a scalar mass squared $m_{h}^{2}$ by
another scalar (fermion: dashed) with the mass $M$. This is a snapshot at
$\mu=2m_{h}$. Its contribution starts from $1-\mu^{2}/m_{h}^{2}$ and quickly
approaches zero.

We can show that all the terms from the same origin are of the same order. For
large $M$, the logarithm is expanded as

$-\log\frac{m_{h}^{2}+M^{2}}{\mu^{2}+M^{2}}=\frac{\mu^{2}-m_{h}^{2}}{M^{2}}-\frac{\mu^{4}-m_{h}^{4}}{2M^{4}}+\frac{\mu^{6}-m_{h}^{6}}{3M^{6}}+\dots.$

Thus the first line in (12) becomes

$\frac{\kappa}{64\pi^{2}}\left[\frac{m_{h}^{4}-\mu^{4}}{2M^{2}}-\frac{m_{h}^{6}-\mu^{6}}{3M^{4}}+\dots\right].$
(13)

Practically, this contribution becomes zero, in the same way the decoupling
theorem applies Appelquist:1974tg . We draw the dependence of the scalar
masses on the Higgs mass squared correction in Fig. 1. Even if a scalar $X$ is
light, its contribution is quickly suppressed for reasonable mass.

Also, for large $m_{t}$, the last line in (12) is

$\frac{y_{t}^{2}}{16\pi^{2}}\sum_{n=1}^{\infty}(-1)^{n}\frac{2n-1}{n+1}\frac{m_{h}^{2n+2}-\mu^{2n+2}}{m_{t}^{2n}}.$
(14)

For considerably large $m_{t}$, the corresponding correction vanishes.222Here,
we mean the case where the external heavy quark mass is given by some other
mechanism than the electroweak Higgs. The top quark cannot be heavy because
its mass is provided by the VEV multiplied by the order one Yukawa coupling.
Again, even for small mass, we see that heavy fermions also easily decouple as
in Fig. 1.

Figure 2: The correction $\delta m_{h}^{2}$ to the Higgs pole mass squared
$m_{h}^{2}$ at renormalization scale $\mu$, by integrating in the top quark
loop and the scalar self-interacting loop (solid). It agrees with the
conventional one-loop correction by top-quark (dotted) Choi:2023cqs . We use
$m_{h}=125.4,m_{t}=173$, all in GeV, and $y_{t}=1,\lambda=0.13$. The
correction $\delta m_{h}^{2}$ is about $-0.39\%$ at 250 GeV.

The Higgs mass correction is insensitive to ultraviolet physics. Only the
high-frequency scalar mass $m_{h}$ is small, and the correction (12) is
sizable.

We plot the correction to the Higgs pole mass squared, as a function of the
scale $\mu$, in Fig. 2. We used the Higgs potential relation for the quartic
coupling $y_{t}v=\sqrt{2}m_{t},m_{h}^{2}=2\lambda v^{2}$.

The one-loop correction to the Higgs mass squared by the top quark is
calculated perturbatively in Ref. Choi:2023cqs . It is also depicted as the
dotted curve in Fig. 2, showing that they match. They are not a priori related
because the present effective field calculation has an additional contribution
by the high-frequency mode of the Higgs scalar, with the new coupling
$\lambda$. It is known that the perturbative one-loop correction by quartic
coupling vanishes Peskin:1995ev .

A similar calculation gives the quartic coupling
$\lambda_{\mu}=\lambda+\delta\lambda(\mu^{2})$ from the reference coupling,

$\begin{split}\delta&\lambda(\mu^{2})=-9\int_{\mu^{2}}^{m_{h}^{2}}\frac{dk^{2}k^{2}}{(4\pi)^{2}}\frac{\lambda_{k}^{2}}{(k^{2}+m_{k}^{2})^{2}}.\\\
&\simeq-\frac{9\lambda^{2}}{16\pi^{2}}\left[\log\frac{2m_{h}^{2}}{\mu^{2}+m_{h}^{2}}+\frac{1}{2}-\frac{m_{h}^{2}}{\mu^{2}+m_{h}^{2}}\right],\end{split}$
(15)

where we approximated $\lambda_{k}$ and $m_{k}$ as before. Note that it is not
necessary to match the coupling $\lambda$ at the Higgs pole mass scale as
before: we can match $\lambda_{k}$ at any scale and run from there, and then
the mass in (15) is the Higgs mass at the matching scale.

Besides the well-known logarithmic running, the last term contains information
on infrared, where light-scalar correction is essential. The change is up to
$0.45\%$ at 250 GeV, justifying the constancy.

We have seen that the natural scale parameter for low energy theory is $\mu$.
By differentiating the total mass with the energy scale $\mu$, we obtain the
renormalization group equation to one-loop order

$\displaystyle\frac{dm_{\mu}^{2}}{d\mu^{2}}$
$\displaystyle=-\frac{3\lambda}{32\pi^{2}}\frac{\mu^{2}}{\mu^{2}+m_{h}^{2}}$
(16)
$\displaystyle\quad-\frac{\kappa}{64\pi^{2}}\frac{\mu^{2}}{\mu^{2}+M^{2}}-\frac{y_{t}^{2}\mu^{2}(\mu^{2}-m_{t}^{2})}{16\pi^{2}(\mu^{2}+m_{t}^{2})^{2}},$
$\displaystyle\mu\frac{d\lambda_{\mu}}{d\mu}$
$\displaystyle=\frac{9\lambda^{2}}{8\pi^{2}}\frac{\mu^{4}}{(\mu^{2}+m_{h}^{2})^{2}}.$
(17)

Using this, we can study multiply coupled equations by various couplings. For
large $m_{h}$ and $m_{t}$, the corresponding fields decouple. For small
$m_{h}$ and $m_{t}$, the right-hand sides become constants independent of
$\mu$, which are commonly used.

Essentially, the same equations can be obtained by differentiating the bare
mass correction (8) and a similar quartic coupling correction with respect to
the upper bound $\Lambda$, instead of the lower bound $\mu$ here. This means
that our exact renormalization group equation has a similar form as that in
Ref. Polchinski:1983gv  (focusing on $\hat{\phi}$)

$\begin{split}\mu\frac{\partial
L}{\partial\mu}&=-\frac{\mu}{2}\frac{\partial}{\partial\mu}\int_{m_{h}<|k|<\mu}\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{k^{2}+m_{k}^{2}}\\\
&\times\left[\frac{\partial L}{\partial\hat{\phi}(k)}\frac{\partial
L}{\partial\hat{\phi}(-k)}+\frac{\partial^{2}L}{\partial\hat{\phi}(k)\partial\hat{\phi}(-k)}\right],\end{split}$
(18)

where $L$ is the momentum-space Lagrangian for the potential $V$ in (4) over
the same regime $m_{h}<|k|<\mu$, including the momentum-conserving delta
function Polchinski:1983gv , using the renormalized couplings. It makes the
partition function

$\displaystyle{\cal Z}$ $\displaystyle=\int{\cal D}\phi{\cal
D}\hat{\phi}\exp\left[-\int_{m_{h}<|k|<\mu}\frac{d^{4}k}{(2\pi)^{4}}\right.$
$\displaystyle\left(\textstyle\frac{1}{2}\hat{\phi}(k)(k^{2}+m_{k}^{2})\hat{\phi}(-k)+J(k)\hat{\phi}(-k)\right)+L\Bigg{]}$
$\displaystyle=\int{\cal D}\phi\exp S_{\text{eff,$\mu$}}$ (19)

invariant under the scale change in $\mu$. It defines the Wilsonian effective
action $S_{\text{eff,$\mu$}}$, containing the mass (10) and the quartic
interaction (15). Conceptually, we track the infrared behavior of the
couplings and see at which energy the correcting fields decouple. Also, using
the energy-dependent renormalized mass gives a more precise result.

In conclusion, the observable parameters of a scalar field are governed by
light fields only. We also calculated the Higgs mass-squared correction from
the pole mass by integrating in its high-energy modes and the top quark.

###### Acknowledgements.

The author is grateful to Jong-Hyun Back and Hans-Peter Nilles for
discussions. This work is partly supported by the grant RS-2023-00277184 of
the National Research Foundation of Korea.

## References

  * (1) K. G. Wilson, Phys. Rev. B 4, 3174-3183 (1971) doi:10.1103/PhysRevB.4.3174; K. G. Wilson, Phys. Rev. B 4, 3184-3205 (1971) doi:10.1103/PhysRevB.4.3184
  * (2) J. Polchinski, Nucl. Phys. B 231, 269-295 (1984) doi:10.1016/0550-3213(84)90287-6
  * (3) K. S. Choi, “On the Observables of Renormalizable Interactions,” [arXiv:2310.00586 [hep-ph]].
  * (4) T. Appelquist and J. Carazzone, Phys. Rev. D 11, 2856 (1975) doi:10.1103/PhysRevD.11.2856
  * (5) M. E. Peskin and D. V. Schroeder, “An Introduction to quantum field theory,” Addison-Wesley, 1995, ISBN 978-0-201-50397-5