# Attention is all you need for boosting graph convolutional neural network

Yinwei Wu The final year project work was carried out under the 3+1+1
Educational Framework at the National University of Singapore (Chongqing)
Research Institute. College of Software Engineering, Sichuan University,
Chengdu, China

###### Abstract

Graph Convolutional Neural Networks (GCNs) possess strong capabilities for
processing graph data in non-grid domains. They can capture the topological
logical structure and node features in graphs and integrate them into nodes’
final representations. GCNs have been extensively studied in various fields,
such as recommendation systems, social networks, and protein molecular
structures. With the increasing application of graph neural networks, research
has focused on improving their performance while compressing their size. In
this work, a plug-in module named Graph Knowledge Enhancement and Distillation
Module (GKEDM) is proposed. GKEDM can enhance node representations and improve
the performance of GCNs by extracting and aggregating graph information via
multi-head attention mechanism. Furthermore, GKEDM can serve as an auxiliary
transferor for knowledge distillation. With a specially designed attention
distillation method, GKEDM can distill the knowledge of large teacher models
into high-performance and compact student models. Experiments on multiple
datasets demonstrate that GKEDM can significantly improve the performance of
various GCNs with minimal overhead. Furthermore, it can efficiently transfer
distilled knowledge from large teacher networks to small student networks via
attention distillation.

## Acknowledgement

This work was carried out under the 3+1+1 Educational Framework at the
National University of Singapore(Chongqing) Research Insitute. The author is
honored to have the support of Prof. Wang Xinchao from the National University
of Singapore, Mr. Wang Xu from Sichuan University and Mr. Jing Yongcheng from
University of Sydney.

###### Contents

  1. 1 Introduction
  2. 2 Related work
    1. 2.1 Graph convolutional neural network
    2. 2.2 Graph Neural network knowledge distillation
    3. 2.3 Attention mechanism
  3. 3 Background
    1. 3.1 Notations
    2. 3.2 Graph nerual network and Graph representation learning
    3. 3.3 Multi-headed self-Attention mechanism
    4. 3.4 Knowledge distillation
  4. 4 Methods
    1. 4.1 Graph knowledge enhancement module
      1. 4.1.1 Reason
      2. 4.1.2 Method
    2. 4.2 Graph knowledge distillation module
      1. 4.2.1 Reason
      2. 4.2.2 Method
  5. 5 Experiments
    1. 5.1 Experimental setup
      1. 5.1.1 Dataset
      2. 5.1.2 Types of GCNs
      3. 5.1.3 Experiment environment
    2. 5.2 The enhancement effect of GKEDM for GCN (RQ1)
    3. 5.3 Analysis of additional parameters GKEDM introduced (RQ2)
    4. 5.4 GKEDM Attention Map Distillation (RQ3)
    5. 5.5 Comparison of GKEDM distillation and related work (RQ4)
  6. 6 Shortcomings and future work
  7. 7 Conclusion

###### List of Figures

  1. 1 The existence of over-smoothing in graph neural networks
  2. 2 GKEDM knowledge enhancement module:GKEDM knowledge enhancement module consists of two phases. In the first phase, GKEDM trains a GCN. In the second stage, GKEDM extracts the trained GCN’s backbone and concats it with the GKEDM knowledge enhancement module for fine-tuning.
  3. 3 Knowledge distillation module of GKEDM:The GKEDM knowledge distillation module will drive the student network to mimic the topology of the teacher network.
  4. 4 Relationship between $\alpha$ setting and distillation effect

###### List of Tables

  1. 1 Notations Table
  2. 2 Dataset statistics
  3. 3 :The effect of GCNII after GKEDM enhancement
  4. 4 :The effect of MoNet after GKEDM enhancement
  5. 5 :The effect of GRAPHSAGE after GKEDM enhancement
  6. 6 :Comparison of number of Param. and Acc.
  7. 7 :Knowledge distillation effect
  8. 8 :Relationship between $\alpha$ setting and distillation effect
  9. 9 :Comparison of GKEDM distillation and related work

##  1 Introduction

Graph Convolutional Neural Networks (GCNs) are powerful representation
learning tools for graph data in non-grid domains. GCNs can learn graph data
structures through generalization of convolution, which is widely used in
computer vision domains. Graph neural networks have made influential
achievements in social networks, recommendation systems, biomedical research
and other domains. Message-passing based GCNs are the prevailing architecture
currently. This type of GCN leverages a message passing algorithm on the graph
to capture node interaction and aggregates node features hierarchically to
extract global features. Through message-passing, nodes obtain information
about themselves and neighbors to update their representations. Message-
passing algorithms typically involve neighbor aggregation, non-linear
activations, and techniques like convolution, pooling, and residual for
complex feature transformations and model optimization. GCNs based on message
propagation mechanism stack multiple graph convolutional layers to achieve a
larger receptive field, enabling better topology capturing of the graph.
However, One major challenge hindering the advancement of GCNs is the issue of
over-smoothing. This refers to a phenomenon in which node representation in a
multi-layer GCN gradually becomes indistinct or overly smoothed as they pass
through layers, making it challenging to distinguish between them. The root
cause of this problem is that in each layer, representations are updated by
averaging or weighted averaging with neighboring nodes’ representations, which
can result in a loss or muddling of information and cause node features to
become increasingly similar. This limitation hinders GCNs from adding extra
layers and therefore restricts their potential to become more expressive
models.

with the deepening of the research, the problem of over-smoothing has
gradually been alleviated [1][2]. As the over-smoothing problem was
alleviated, researchers began gradually increasing the number of parameters of
the GCNs in order to boost their expressive capabilities, just as in the image
and text domains. This trend makes it a challenge to deploy graph neural
networks under resource and time constraints.

Knowledge distillation was proposed by [3] in 2015, and has now become the
mainstream method of model compression. Knowledge distillation provides
additional supervision signals for the training of the student network to
improve its performance. The additional supervision signals come from a
trained teacher model with a large number of parameters, which extracted the
hidden distribution of the complex training data. Initially, knowledge
distillation relied solely on soft target probabilities imparted by teachers
to facilitate the learning of student networks. Recently, to further
investigate the potential of knowledge distillation, many studies are no
longer confined to the soft target probabilities. Instead, they also utilize
the intermediate layers of the model for distillation [4][5].

Knowledge distillation in traditional neural networks for image and text data
has been extensively studied. However, research on distilling knowledge in
graph neural networks is still in its infancy [6][7][8]. Graph data contains
unique information, such as its topological structure, and is distinct from
grid data. To improve distillation efficiency, careful selection and design of
the information to be distilled is necessary.

In this work, a novel plug-in module for GCNs named GKEDM (Graph Knowledge
Enchance and Distillation Module), is proposed. GKEDM is capable of replacing
the output layer of GCN. Its knowledge enhancement module employs an attention
mechanism for local topology information aggregation to boost node
representation, thus improving network performance.

Simultaneously, to alleviate the computational burden resulting from
parameters, a custom attention graph knowledge distillation method that
suitable for GKEDM modules is introduced. This method enables the teacher
network to transfer local topology to the student model, achieving the goal of
distillation.

To demonstrate the universality of GKEDM, I performed enhancement and
knowledge distillation on different kinds of GCN. Experimental results across
multiple datasets demonstrate that GKEDM significantly enhances the
capabilities of graph neural networks and optimizes the efficiency of
knowledge distillation.

The contributions of this paper are:

  * •

A plug-in GCN enhancement module based on graph local topology attention is
proposed, named GKEDM. It can be directly applied to GCNs and improve its
performance.

  * •

A novel knowledge distillation method suitable for GKEDM is introduced, which
can effectively transfer the topology knowledge of GCN from teacher network to
student network.

  * •

The experimental results show that the GKEDM is versatile and effective in
various GCNs and datasets.

##  2 Related work

This work is related to graph convolutional neural networks, knowledge
distillation on graph convolutional neural networks, and attention mechanisms.

###  2.1 Graph convolutional neural network

In recent years, the domain of graph neural networks(GNN) has gained
significant research attention owing to its powerful processing capabilities
on non-grid data. After the convolution operation on the graph was defined
[9], an increasing number of convolution-based approaches are being used to
enhance the representational capacity of GNNs. That specific type of GNN is
referred to as GCN. Graph Attention Network(GAT) [10] is an attention
mechanism-based GCN designed for processing graph data. Compared to other
GCNs, GAT utilizes a trainable inter-node attention mechanism to evaluate the
relevance of each node to the target node. It then performs weighted
aggregation of representations of neighboring nodes to better capture the
local topology. his attention mechanism adapts to the representations of
neighboring nodes, making GAT more flexible and comprehensible while
processing graph data. GraphSAGE [11] leverages self-supervised learning to
update a node’s embedding based on its neighbors, providing flexibility in
handling diverse graph data types.With neighbor sampling, it can efficiently
process large-scale graph data. GCNII [12] adopts residule block to prevent
over-smoothing while training, thereby enabling the stacking of more GCN
layers. MONET [13] incorporates a learnable distance-based kernel function
that transforms similarity between nodes into Gaussian Kernel Functions at
multiple scales. This approach retains information at various levels,
resulting in improved node representations. Rather than designing a new
structure of GCN, this work focuses on improving the performance of the
existing GCNs, and performing knowledge distillation and model compression on
them.

###  2.2 Graph Neural network knowledge distillation

The data to be processed by convolutional neural networks(CNN) performing on
grid domain data is often a fixed-size matrix, such as image data. Different
from traditional CNN, the GCNs need to deal with the non-regular graph data
composed of nodes and edges. With additional features that grid domain data
does not have, such as node features, graph topology, communities, etc., GCNs
usually need to be carefully designed. Taking into account the aforementioned
dissimilarities, the knowledge distillation approach for GCNs must also
address the distillation of supplementary graph data structure captured. GCN
knowledge distillation was first proposed in [6]. Yang et al. captured the
realationship between nodes by introducing an LSP(Local Structure preserving)
module to achieve local topology distillation in graphs. More specifically,
LSP uses the distance of the node representation in the kernel space to
measure the similarity between nodes, and regards the similarity of the first-
order neighborhood of nodes as a distribution. By letting the student network
mimic the node neighborhood distribution of the teacher network, the purpose
of distillation topology knowledge can be achieved. Subsequently, other
methods of graph neural network knowledge distillation were proposed. In [8],
He et al. used an adversarial way to perform model compression, named
GraphAKD. In GraphAKD, the student network as a generator needs to produce
output similar to the teacher network to fool the discriminator, and the
discriminator needs to distinguish between the output of the student and the
output of the teacher. Contrastive learning distillation has been proposed in
the image domain [14], and contrastive learning in the graph neural network
domain has also been used as a way of representation learning [15] to prove
effective. CKJ et al. [7] combined graph neural network distillation and
contrastive learning to propose G-CRD. Through contrastive learning, G-CRD
makes the representation of the student network node approximate the
representation of the corresponding node of the teacher network, and at the
same time, it is far away from the representation of other nodes of the
teacher network, implicitly preserving the relationship between nodes, and
achieving the purpose of topology distillation. The knowledge distillation
approach introduced in this research differs from previous approaches, as it
focuses on a specific architecture for graph neural networks with GKEDM
enhancement.

###  2.3 Attention mechanism

Attention mechanism (AM) is an important technique for enhancing the
expressive ability of neutral network models. Attention allows for varying
weights to be assigned to different segments of input data. This enhances
neural networks’ ability to manage complex data structures and multiple
sources of information. The AM is actually reflected in the very early gate
mechanism, but it received extensive attention after the Transformer [16]
framework came out. Compared to traditional convolution, AM has stronger
relationship capturing abilities and significant potential applications in the
era of big data. Transformer was first applied to the natural language
processing domain, and then migrated to the image domain [17], which has also
proved to be very effective. All the current giant models with amazing results
are basically evolved based on Transformer architecture.

In the graph neural network domain, the attention mechanism was first used in
GAT [10]. Unlike Transformer, the GAT does not use the framework of Query,
Key, and Value, but the idea based on attention still achieves good results.
Since the GAT is implemented based on the message passing mechanism, there can
be problems with over-smoothing when training large deep models. In order to
solve the oversmoothing problem and design a Transformer suitable for graph
data, Ying et al. abandoned the architecture based on the message passing
mechanism, and carefully designed a very important feature in the Transformer
in the proposed Graphormer [2], which is positional encoding. Positional
encoding has been proposed in the Transformer of the original version for
natural language processing. Since Transformer does not have explicit position
information, and order information is crucial in text sequences, a mechanism
is needed to represent the position of each input. To solve this problem,
Transformer added positional encoding to explicitly identify the position
information of the input for the model. Positional encoding has been proved to
have a significant influence on the results in the design of Transformer, so
subsequent Transformer have also adopted a variety of positional encodings
suitable for different data structures. Graphormer can excel in graph
structures by creating appropriate positional encodings to denote the roles,
locations, and other characteristics of nodes within the graph. However,
Graphormer, as a derivative architecture of Transformer, also has the most
important problem, that is computational complexity. And because of the
particularity of the graph structure, this problem becomes more serious.Each
node must query and compute with all nodes in the graph, resulting in a
computational complexity of $O(n^{2})$. This challenge presents a major
obstacle for the application of Graphormer to graph datasets containing tens
of thousands of nodes.

This work employs an attention mechanism that leverages the strengths of GAT
and Graphormer and incorporates a self-attention mechanism based on local
structure, thereby reducing computational load. This work incorporates a self-
attention mechanism utilizing Query, Key, and Value architectures, while also
implementing positional encoding based on the graph Laplacian matrix and
random walk to enhance the model’s expressive capability.

##  3 Background

The purpose of this work is to enhance and distill GCNs based on attention
mechanism. This section first introduces some basic principles and algorithms
of graph representation learning based on message passing mechanism, and is
followed by an introduction related to the self-attention mechanism and
knowledge distillation.

###  3.1 Notations

Due to the large number of symbols involved, this subsection provides a
comprehensive listing of all mathematical symbolic representations employed in
this work, as well as their corresponding descriptions, for the sake of
clarity and convenience. The notations are shown in Table 1:

Notation | Description
---|---
$G=(\mathcal{V},E)$ | A graph $G$ containing the set of edges $E$ and the set of nodes $\mathcal{V}$
$v_{n}\in\mathcal{V}$ | Denotes the node $n$
$e_{ij}\in E$ | Denotes the edge connecting nodes $i$ to $j$
$H_{i}^{T},H_{i}^{S}$ | Denotes all the node representations at the output of layer $i$ of the teacher GCN and
all the node representations at the output of layer $i$ of the student,
respectively
$h_{n}^{i}$ | Denotes the representation of node $n$ at layer $i$
$z_{n}$ | Denotes the logits of node $n$ predicted by the GCN
$y_{n}$ | Denotes the label of node $n$
${y_{n}^{pred}}$ | Denotes the prediction result of the network for node $n$
$x_{n}$ | Denotes the initial features of node $n$
$N(v_{n})$ | Denotes 1-hop neighbors of node $n$
$Q,K,V$ | Denotes $Query,Key,Value$ in attention mechanism
$d$ | Denotes the length of the node embedding vector
Table 1: Notations Table

###  3.2 Graph nerual network and Graph representation learning

Graph neural networks (GNNs) are neural networks that can generate node
embeddings based on node and edge features as well as graph topology. These
embeddings can be further used for various downstream tasks. Types of tasks
that GNN can perform include node prediction, graph prediction, edge
prediction, and relationship prediction. This work primarily focuses on node
classification of graph convolutional neural networks (GCN). Therefore, this
section delves into the GCN and its message passing mechanism, which is a
fundamental algorithm for performing classification predictions.

GCN receives graph that composed of edge set $E$ and node set $\mathcal{V}$ as
input. For each node $v_{n}\in\mathcal{V}$ there exists an initial
$d-dimensional$ input feature vector $x_{n}$, and in general, this initial
feature is the input to layer 0 of the GCN. Also, each node has a
corresponding label $y_{n}$, and 1-hop neighbor node set $N(v_{n})$. For the
node classification task, the GCN needs to predict a label ${y_{n}^{pred}}$
for each node based on its feature information and the topology of the graph,
and make ${y_{n}^{pred}}$ and $y_{n}$ as similar as possible in terms of
metrics. A message passing GCN consists of a stack of convolutional layers,
each of which can be viewed as an encoder. In each layer of the GCN, the nodes
embed the topology in their own information by collecting information from
their neighbors. This stage is known as $AGGREGATE$. Then the node enters the
second $UPDATE$ phase to generate and update its own node representation in
the next layer. The above aggregation-update process can be expressed in
mathematical language as the equation (1):

$\displaystyle message_{N(v)}^{i}=AGGREGATE^{(i)}({h_{u}^{i}|\forall u\in
N(v)})$ (1) $\displaystyle
h_{v}^{i+1}=UPDATE^{(i)}(h_{v}{{}^{i}},message_{N(v)}^{i})$

In the above equation, $AGGREGATE$ is an aggregation function with
permutation-invariant, $message$ is the information of the local neighborhood
that the node aggregated, and $h_{v}^{i+1}$ is the input representation of the
node in the next layer, specifically, $h_{v}^{0}=x_{v}$. After a total of $i$
layers of message passing convolutional layers are stacked, the node can
aggregate the information of the nodes in the $i$-order neighborhood and thus
has a more powerful representation. For the node classification task, the GCN
ultimately passes the node’s final representation to a linear classifier, such
as a multilayer perceptron (MLP). This step calculates the probability
distribution $z_{v}$ for each class and determines the loss $L_{CE}$ using the
cross-entropy loss function as stated in Eq. (2):

$\displaystyle z_{v}$ $\displaystyle=MLP(h_{v}^{i})$ (2) $\displaystyle
L_{CE}$ $\displaystyle=H(y_{v},z_{v})$

$H$ denotes the cross-entropy loss function. After the training is completed
by gradient descent, the prediction results of the nodes can be obtained by
$softmax(z_{v})$.

###  3.3 Multi-headed self-Attention mechanism

Self-Attention mechanism based on $Key,Query,Value$ has received extensive
attention after being proposed in Transformer[16]. Its main principle is to
map the set of input feature vectors $H\in\mathbb{R}^{n\times d}$ to three new
sets of features located in different expression spaces ,
$K\in\mathbb{R}^{n\times d},Q\in\mathbb{R}^{n\times d},V\in\mathbb{R}^{n\times
d}$, by a function as shown in Eq. (3):

$\displaystyle K=f_{k}(H),$ (3) $\displaystyle Q=f_{q}(H),$ $\displaystyle
V=f_{v}(H),$

For any two input nodes, the similarity is computed by point multiplication
using the $Query$ features of the source node and the $Key$ features of the
target node. Once the similarity between the target node and all the source
nodes is obtained, the target node weights and sums the value features of
these source nodes based on the calculated similarity to produce its own new
representation. The mathematical expression is shown in the Eq. (4):

$\displaystyle Attention(Q,K,V)=Softmax(QK^{T}/\sqrt{d})V,$ (4)

$QK^{T}$ is a $n\times n$ matrix that represents the similarity between two
two nodes, and $\sqrt{d}$ is a normalization factor. The final resulting
$Attention(Q,K,V)$ matrix has the same shape as $V$, so that the Transformer
can improve the expressiveness of the model by stacking multiple identical
attention layers.

In order to enable the model to observe more representation space, Multi-head
Self Attention (MSA) is one of the main extensions to enhance the
representational power of Transformer models. The MSA is based on the idea of
splitting, where the input vector is partitioned into multiple subvectors and
each subvector is assigned a separate attention head. Each attention header
contains an independent set of attention weight matrices for computing the
self-attention. Finally, the self-attention results of all heads are
aggregated to form the final output vector. For a MSA using $k$ heads, the
model will use $K$ sets of different functions to obtain $K$ sets $Q,K,V$ and
perform the self-attentive mechanism separately. Each head will provide a
matrix in the shape of $V$. Finally, for all the results of self-attention, a
linear mapping or averaging can be used to obtain the output. This approach
can enhance the performance of the Transformer to some extent.

###  3.4 Knowledge distillation

Knowledge distillation is a technique for enhancing the performance of student
models through the transfer of advanced and complex knowledge from large and
high-performance teacher models. By feeding the output of the teacher model as
a supervised signal to the student model, the student model can leverage the
experience and knowledge of the teacher model to enhance its ability to learn
the distinguishing features and patterns of the task, ultimately improving
overall performance. The knowledge distillation first proposed in [3] by
having the student output label probabilities to match the teacher output soft
label probabilities, with a specific loss function formula of Eq. (5):

$\displaystyle L_{KD}=H(softmax(z^{T}/T_{1}),softmax(z^{S}/T_{1}))$ (5)

where $H$ denotes the cross-entropy loss or KL divergence, and $T_{1}$ denotes
the distillation temperature, which allows for smoother labels for
distillation. $z^{T}$ and $z^{S}$ are the labeling probabilities of teachers
and students, respectively.

##  4 Methods

In this section, we will briefly introduce the motivation for proposing GKEDM
and present the principles and specific algorithms of GKEDM.

###  4.1 Graph knowledge enhancement module

####  4.1.1 Reason

Message passing based GCN embeds topology into the representation vector of
nodes by aggregation of multiple convolutional layers. However, overstacking
the number of convolutional layers can cause over-smoothing, i.e., the
representation of all nodes tends to be the same and cannot be differentiated.
One of the reasons for this phenomenon is that GCNs indiscriminately aggregate
information from node neighbors, and some information from important neighbors
may be underestimated, while some noisy information from completely unrelated
neighbors may be amplified. Over-smoothing phenomenon becomes more and more
significant as the number of layers of the GCN increases, which in turn leads
to degradation of the model performance.

From the figure (1), it can be seen that in a variety of GCNs implemented
based on the message passing mechanism, the prediction performance reaches
saturation and declines as the number of layers increases.Therefore, it is
more common to build shallow graph GCNs with only 2-3 layers. While specially
designed graph GCNs, such as GCNII, address the over-smoothing problem by
utilizing a residual connection structure similar to that in image domains,
the indiscriminate aggregation of neighboring node information still adversely
impacts model performance to a certain extent. The additional computational
effort required by this approach cannot be ignored. It is crucial to find ways
to improve the performance of GCNs while controlling the amount of parameters.

Figure 1: The existence of over-smoothing in graph neural networks

The MSA, which is widely used in Transformer, can calculate the similarity of
two nodes’ representations, and then selectively choose the information passed
by neighbors for weighted reception. By selectively receiving different
information through the attention mechanism, nodes are able to filter out
useful information and discard noisy information, which is more conducive to
the functional differentiation of nodes, i.e., making different types of nodes
more distinguishable. It is beneficial for downstream node classification
tasks.

However, given the unique characteristics of graph data’s non-grid structure,
applying self-attention mechanisms to graph neural networks requires special
consideration. If applies the attention mechanism to all node pairs and
aggregates the global topology, although it can increase the expressiveness of
the model, the $O(n^{2})$ time complexity will be a huge burden. Therefore,
GKEDM chooses to trade-off the expressiveness and computational overhead by
choosing to aggregate node information in the first-order neighborhood of
nodes with MSA. Unlike GAT, GKEDM adopts the paradigm of Transformer multi-
headed attention mechanism and adds location encoding to enhance the
expressiveness and differentiation of the nodes.

####  4.1.2 Method

The core goal of GKEDM is to enhance and coalesce node attribute information
and graph topology information on the backbone of various types of GCN
networks that have been trained to obtain a better node representation. the
algorithm flow of the knowledge enhancement part of GKEDM is shown in Figure
3.

Figure 2: GKEDM knowledge enhancement module:GKEDM knowledge enhancement
module consists of two phases. In the first phase, GKEDM trains a GCN. In the
second stage, GKEDM extracts the trained GCN’s backbone and concats it with
the GKEDM knowledge enhancement module for fine-tuning.

GKEDM’s knowledge enhancement module is a two-stage algorithm. In the first
stage, GKEDM first pre-trains a GCN, as shown in the upper part of Figure 3.
We assume that the lower layers of GCN are mainly responsible for learning the
representation of nodes, while at higher layers it is more concerned with
aggregating the topology into the graph node representation. Thus the main
task at this stage is to make the lower convolutional layers of the GCN learn
the basic representation of a node through optimization.

After learning the basic representation of nodes, GKEDM replaces the last
convolutional layer of the GCN with a MSA based layer in order to allow the
network to better aggregate the information of neighboring nodes in the higher
convolutional layers and thus enhance the capture of topology. The MSA based
graph convolutional layer can perform more effective message aggregation.
Instead of choosing to add an additional MSA based layer after the last layer
of the convolutional layer of the trained GCN, GKEDM chose to perform a
replacement, an approach with two considerations:

  * •

The excessive stacking of graph convolutional layers may cause information
loss and smoothing of node features. Therefore, reducing the number of layers
of graph convolutional neural network can avoid the over-smoothing phenomenon.

  * •

Reducing the number of layers of graph convolution layers can reduce the
computational burden to some extent.

To apply the GKEDM module on the GCN, for a pre-trained GCN with $k+1$ layers,
a node feature matrix $H_{k}\in\mathbb{R}^{n\times d}$ is generated at the
$k$th layer after feeding graph $G$. At this point $H_{k}$ can be considered
to have learned the basic representation of each node. To further
differentiate the nodes and enhance their distinguishability, thus
facilitating the downstream node classification task, GKEDM computes a
positional encoding (PE) for each node. PE has proven to be very effective for
applications in Transformer in the image and text domains. The principle of
its action is to allow the network to distinguish the positional information
of the input. In graph, nodes may have different responsibilities,
communities, etc., and PE can provide the network with effective priori by
using this information. GKEDM uses the Laplace position code [18]. The
Laplacian matrix is a very important mathematical tool to measure the graph
properties and it can represent the location characteristics of the nodes in
the graph. Laplacian PE is done by extracting the top $m$ smallest non-
repeating eigenvalues in the graph Laplacian matrix as the eigenvectors for
location coding. After combining the above position encoding, the new input
$\hat{H_{k}}$ of GKEDM is computed as Eq. (6):

$\displaystyle\hat{H_{k}}=H_{k}+f(PE(G))$ (6)

$PE$ is the Laplace PE, and $f$ is a linear mapping layer that maps the PE to
the $d$-dimension node hidden vector. After obtaining $\hat{H_{k}}$, the MSA
layer is computed using Eq. (3) and Eq. (4) to get $Attention$, and residual
connection is used to stabilize the performance.

$\displaystyle K$ $\displaystyle=f_{k}(\hat{H_{k}}),$ (7) $\displaystyle Q$
$\displaystyle=f_{q}(\hat{H_{k}}),$ $\displaystyle V$
$\displaystyle=f_{v}(\hat{H_{k}}),$ $\displaystyle A$
$\displaystyle=Softmax(QK^{T}/\sqrt{d})$ $\displaystyle H_{k+1}$
$\displaystyle=AV+\hat{H_{k}},$

Finally, the GKEDM can be optimized by the loss function of Eq. (2) and the
stochastic gradient descent algorithm.

###  4.2 Graph knowledge distillation module

####  4.2.1 Reason

Although the knowledge enhancement module of GKEDM is able to improve the
performance of the GCN through MSA, it introduces additional parameters
proportional to the length $d$ of the hidden vector of the node
representations due to the computation of $Q,K,V$. When $d$ is large, adding
GKEDM for knowledge enhancement introduces a large computational burden.

Knowledge distillation is a common model compression technique applied in
neural networks. Since the classification performance of GCN has been greatly
improved by GKEDM, knowledge distillation can be performed at the cost of a
small performance loss, and a trade-off between computational speed and
accuracy can be made to achieve both improved accuracy and reduced model size.

####  4.2.2 Method

There’s some research on how to apply knowledge distillation to graph neural
networks [6][8][7]. In this work, a customized distillation method for GKEDM,
named _Attention Map Distillation (AMD_), is used, which is different from the
above method. AMD is closer in principle to [6]’s Local Structure Perserving
(LSP) structure, so we first give a brief introduction to the LSP structure.

The LSP structure treats the 1-hop neighborhood of a node as distribution, and
the purpose of its distillation is to train the student network to mimic the
1-hop neighborhood distribution learned by the teacher network. To be more
precise, as GCN goes deeper and incorporates topological information into node
representations, the distribution of node neighborhoods can be illustrated as
the similarity between neighboring node representations and central node
representations in the hidden space. In the LSP, the similarity between nodes
is calculated with three kernel functions:

$\displaystyle K(h_{n}^{i},h_{n}^{j})=\left\\{\begin{aligned}
&((h_{n}^{i})^{\top}h_{n}^{j})^{d},&Poly\\\
&e^{-\frac{1}{2\sigma^{2}}}||h_{n}^{i}-h_{n}^{j}||^{2}&RBF\\\
&(h_{n}^{i})^{\top}h_{n}^{j},&Linear\end{aligned}\right.$ (8)

The kernel function can map the low-dimensional features into the high-
dimensional feature space, thus making the samples that can not be separated
in the low-dimensional space linearly separable in the high-dimensional space.
Among the above three kernel functions, the RBF kernel function is the most
widely used, which can map the samples into an infinite dimensional feature
space with good nonlinear fitting ability and adaptability. LSP uses the KL
divergence to measure the difference in the 1-hop neighborhood distribution
around each node for teachers and students, and uses this as a loss function
for optimization:

$\displaystyle L_{LSP}=\sum_{i\in\mathcal{V}}D_{KL}(\mathop{softmax}_{i,j\in
E}(K((h_{n}^{i})^{S},(h_{n}^{j})^{S}))||\mathop{softmax}_{i,j\in
E}(K((h_{n}^{i})^{T},(h_{n}^{j})^{T})))$ (9)

It can be seen that the LSP uses a predefined kernel function to measure the
similarity of two nodes. It may have some limitations:

  * •

The assumption of using kernel functions is that the representations of two
similar nodes in the kernel function space are also similar, which is not
guaranteed. Similar node representations do not necessarily remain similar
after the mapping of kernel functions. The kernel function approach lacks
flexibility.

  * •

The teacher network and the student network may have difference in parametric
size, so the student network may limit its own exploration by completely
imitating the structure of the teacher network in the kernel function space.

To this end, in comparison to the LSP structure that uses a predefined
unlearnable kernel function to characterize the similarity between
representations, the AMD proposed in this paper uses learnable parameters to
measure the distance between nodes. AMD improves the effectiveness of
knowledge distillation by training to fit a more suitable measure of node
similarity. A similar approach to attention graph distillation was proposed in
[19], but its role is in the field of natural language processing, while this
paper is the first work to apply attention distillation to graph neural
network distillation.

The knowledge distillation module of GKEDM is shown in the figure:

Figure 3: Knowledge distillation module of GKEDM:The GKEDM knowledge
distillation module will drive the student network to mimic the topology of
the teacher network.

The GCN for knowledge enhancement using GKEDM contains a MSA layers introduced
by GKEDM. and the attention map generated by its MSA layer can be used as
additional supervised information for knowledge distillation on the student
network after the teacher network has completed training. The attention maps
$A^{T}$ and $A^{S}$ of teacher and student respectively can be obtained by the
equation (6), which are representations of the similarity of nodes to their
neighbors, and thus a description of the node neighborhood structure. By
driving the student network to mimic the neighborhood structure of the teacher
network, GKEDM is able to achieve the purpose of attention distillation.
Different from LSP, the attention map of this method is obtained by $Query$,
$key$ computation containing learnable parameters, and gains more flexibility
through the MSA. This approach can overcome the distillation instability
problem of teacher and student networks due to size mismatch. In particular,
we train the teacher and student attention graphs by minimizing the KL
divergence of:

$\displaystyle L_{A}=L_{KL}(A^{T}||A^{S})$ (10)

With attention graph distillation, $Key$ and $Query$ pair information between
teacher network nodes is transferred to students. However, in addition to Key
and Query information, $Value$ are also very important representations of node
descriptions in the MSA. In order to reach a more efficient distillation,
Value-Value relations between teacher node pairs are also delivered as an
additional attentional information for students to learn:

$\displaystyle L_{VR}=L_{KL}((V^{T})^{\top}V^{T}||(V^{S})^{\top}V^{S})$ (11)

$V^{T}$ and $V^{S}$ are derived from the teacher and student, respectively.
The above equation measures the dissimilarity between node pairs by the
$Value-Value$ pairs. It is worth noting that this approach can be applied to
Key-Key pairs and Query-Query pairs, which will be compared in the
experimental section of this paper. The total loss function of the attention
distillation is finally expressed as:

$\displaystyle L_{attention}=L_{A}+L_{VR}$ (12)

And the final loss function of the student network training is the cross-
entropy loss plus the attention loss, which is expressed by the Eq. 13

$\displaystyle L=L_{CE}+\alpha L_{attention}$ (13)

$\alpha$ is the weight of the attention loss function, which indicates the
intensity of attention distillation.

##  5 Experiments

GKEDM is proposed to enhance the GCN for node knowledge extraction and to
compress GCN’s knowledge to a smaller network by distillation. In order to
verify the feasibility and generality of GKEDM, we set up the experiments to
answer the following questions.
RQ1:Can GKEDM work for different types of GCN on datasets of different sizes
and tasks?
RQ2:Does GKEDM’s improvement in performance come at the cost of additional
parameters?
RQ3:Does GKEDM’s attention map distillation really work?
RQ4:How does GKEDM’s attention map distillation compare to other distillation
methods?
In this section, we first introduce the dataset, the GCN and the experimental
environment used, followed by experiments to answer each of the above four
questions.

###  5.1 Experimental setup

####  5.1.1 Dataset

We conducted experiments on a series of datasets containing different numbers
of graphs and nodes for the node classification task. These datasets include
both multi-class tasks and multi-label tasks. The statistical information of
the datasets is summarized in the table (2). More detailed information is as
follows .

  * •

PPI(Protein-Protein Interaction) [20] is a protein interaction dataset. Each
protein is represented as a node and the interactions are represented as
edges. Its contains a total of 24 graphs and is a multi-label classification
task.

  * •

The task of the FLICKR [21] dataset is to classify images on the web by their
description and properties.

  * •

CORA FULL [22] is an expanded CORA dataset. It contains a graph representing
paper citation relationships, with each node representing a paper. The task of
this dataset is to predict the kind of papers.

Dataset | Number of nodes | Number of node’s feature | Task
---|---|---|---
PPI | 2,372 | 50 | Multi-label(121 labels)
FLICKR | 89,250 | 500 | Multi-class(7 classes)
CORA FULL | 19,793 | 8710 | Multi-class(70 classes)
Table 2: Dataset statistics

####  5.1.2 Types of GCNs

In addition to experiments on different datasets, to demonstrate that GKEDM
does not work only for specific GCN, we also conducted experiments on
different GCNs:

  * •

GCN [9] is the most classical semi-supervised graph neural network model. It
learns node representations by defining convolutional operations on graph data
structures.

  * •

GCNII [12] avoids the over-smoothing problem existing in GCNs by residual
connections, and thus can be stacked to higher layers.

  * •

GRAPHSAGE [11] is a GCN that learns node representations by sampling and
aggregating the neighbors of nodes.

  * •

MONET [13] is a network that can be used on non-Eulerian domains, and is
therefore suitable for the graph. It enhances graph networks by using kernel
functions.

####  5.1.3 Experiment environment

All experiments were performed on a LINUX system using pytorch framework
version 1.12.1+cu116. The GNN framework is DGL, version 1.0.2+cu116, and the
hardware used is a RTX3090.

###  5.2 The enhancement effect of GKEDM for GCN (RQ1)

To improve the performance of GCNs on node classification tasks, this work
introduces an enhancement module, GKEDM. GKEDM aims at weighting the
aggregated node neighborhood information and updating the basic representation
of nodes by introducing an attention mechanism. However, datasets of different
sizes and tasks may differ in the structure of the graph, and the semantics of
the node features. Also, since the principles of different kinds of GCNs are
different, the node representations they learn may also differ. To verify
whether the attention mechanism can ignore these differences and achieve its
original purpose of aggregating valid information, this experiment uses
different GCNs trained on different datasets and compares the classification
accuracy of the network without GKEDM with that of the network with GKEDM.
Since the PPI dataset is a multi-label dataset, the metric used is the
F1-score. The experimental results are shown in the table (3)(4)(5):

Model | Dataset | Layers | Params(M) | Original Acc. | Acc. with GKEDM | Impv.
---|---|---|---|---|---|---
GCNII | PPI | 16 | 2.14 | 0.82 | 0.97 | 0.15
GCNII | PPI | 16 | 0.54 | 0.55 | 0.85 | 0.30
GCNII | CORA FULL | 1 | 0.28 | 0.64 | 0.85 | 0.21
GCNII | CORA FULL | 8 | 3.3 | 0.70 | 0.89 | 0.19
Table 3: :The effect of GCNII after GKEDM enhancement Model | Dataset | Layers | Params(M) | Original Acc. | Acc. with GKEDM | Impv.
---|---|---|---|---|---|---
MoNet | PPI | 3 | 0.03 | 0.80 | 0.93 | 0.13
MoNet | PPI | 3 | 0.54 | 0.82 | 0.97 | 0.15
MoNet | CORA FULL | 3 | 5.7 | 0.65 | 0.86 | 0.21
MoNet | CORA FULL | 2 | 0.28 | 0.66 | 0.81 | 0.15
Table 4: :The effect of MoNet after GKEDM enhancement Model | Dataset | Layers | Params(M) | Original Acc. | Acc. with GKEDM | Impv.
---|---|---|---|---|---|---
GRAPHSAGE | PPI | 5 | 0.24 | 0.76 | 0.92 | 0.16
GRAPHSAGE | PPI | 3 | 0.03 | 0.66 | 0.85 | 0.19
GRAPHSAGE | CORA FULL | 2 | 0.14 | 0.69 | 0.84 | 0.15
GRAPHSAGE | CORA FULL | 2 | 0.28 | 0.65 | 0.86 | 0.21
Table 5: :The effect of GRAPHSAGE after GKEDM enhancement

From the experimental results, it can be seen that the performance of the GCNs
on the node classification tasks is significantly enhanced after applied
GKEDM, and the accuracy improvement of up to 30% can be achieved. The
enhancement effect is demonstrated across various GCNs and datasets,
highlighting the generality of GKEDM for GCN enhancement. It is worth noting
that, as shown in Table 3 for experiments on GCNII, different degrees of
performance enhancement are observed even for same models with large
differences in the number of layers and parametric quantities. It indicating
that the GKEDM is a general module that works effectively for features of
different scales and semantics.

###  5.3 Analysis of additional parameters GKEDM introduced (RQ2)

For neural networks, increasing the number of parameters can improve the
fitting ability of the model, which is one of the most direct ways to improve
the performance. The additional trainable parameters allow the model to
improve the learning and to capture the information that could not be captured
originally. The knowledge enhancement module of GKEDM uses an MSA layer to
replace the last layer of the original GCN. The attention layer maps node
representations to $Query$, $Key$, $Value$, thus introducing additional
parameters. Although the GKEDM-enhanced model has better performance, the
introduction of additional parameters makes the number of parameters more
significant compared to the original model. Therefore, it is not convincing to
compare only the accuracy of the original model with that of the GKEDM
enhanced model. The purpose of this experiment is to demonstrate that the
improvement of GKEDM for model performance does not rely on the introduction
of additional parameters, but on more efficient parameter utilization. This
experiment uses the number of parameters and accuracy of the model as metrics
for comparison, and the experimental results are shown in Table 6:

Model | Dataset | Total Param.(M) | With GKEDM | Acc.
---|---|---|---|---
GCNII | PPI | 2.14 | No | 0.82
GCNII | PPI | 0.66 | Yes | 0.85
GCNII | CORA FULL | 3.3 | No | 0.64
GCNII | CORA FULL | 0.28 | Yes | 0.85
GRAPHSAGE | PPI | 0.24 | No | 0.76
GRAPHSAGE | PPI | 0.14 | Yes | 0.85
GRAPHSAGE | CORA FULL | 0.29 | No | 0.65
GRAPHSAGE | CORA FULL | 0.14 | Yes | 0.84
MoNet | PPI | 0.54 | No | 0.82
MoNet | PPI | 0.22 | Yes | 0.93
MoNet | CORA FULL | 5.7 | No | 0.65
MoNet | CORA FULL | 0.28 | Yes | 0.84
Table 6: :Comparison of number of Param. and Acc.

As illustrated in Table 6, we performed a comparison of performance and number
of parameters for the same GCN with or without the addition of GKEDM. From the
results, it can be seen that the GCN with GKEDM is able to surpass the
performance of the original GCN with a small number of parameters, which
proves that GKEDM does not rely on number of parameters to improve the
expressive power of the model. From the comparison of MoNet on the dataset
CORA FULL below Table 6, even a small model with only one twentieth of the
number of parameters can outperform a large model in the final performance
through the enhancement of GKEDM. This experiment illustrates the ability of
GKEDM to improve the parameter utilization of GCNs and the effectiveness of
the GKEDM for GCNs enhancement.

###  5.4 GKEDM Attention Map Distillation (RQ3)

The main purpose of this experiment is to verify the effectiveness of GKEDM
attention map distillation. The task of knowledge distillation is mainly to
improve the performance of the small student model by passing the dark
knowledge of the large teacher model as an additional supervisory signal.
Therefore, the number of teacher network participants should be more than the
student network and the performance should be appropriately separated from the
student network. Therefore, this experiment was chosen to be conducted on
FLICKER and PPI where the distillation effect could be seen. Meanwhile, in
order to try the effect of different combinations of attentional collation, we
also tried Key-Relation and Query-Relation in addition to distillation of
Value-Relation. The experimental results are shown in the following figure:

Model | Dataset |  Teacher Param.(M) | Teacher Acc. |  Student Param.(M) |  Student original Acc. | Distill Acc. | Type
---|---|---|---|---|---|---|---
GraphSAGE | PPI | 0.74 | 0.9452 | 0.15 | 0.8624 | 0.8755 | a+v
GraphSAGE | PPI | 0.74 | 0.9452 | 0.15 | 0.8624 | 0.8750 | a+v+q+k
GraphSAGE | FLICKR | 0.4 | 0.8183 | 0.2 | 0.6983 | 0.7196 | a+v
GraphSAGE | FLICKR | 0.4 | 0.8183 | 0.2 | 0.6983 | 0.7147 | a+v+q+k
MoNet | PPI | 1.04 | 0.9753 | 0.33 | 0.9264 | 0.9399 | a+v
MoNet | PPI | 1.04 | 0.9753 | 0.33 | 0.9264 | 0.9399 | a+v+q+k
MoNet | PPI | 1.04 | 0.9753 | 0.21 | 0.9023 | 0.9281 | a+v
MoNet | PPI | 1.04 | 0.9753 | 0.21 | 0.9023 | 0.9302 | a+v+q+k
Table 7: :Knowledge distillation effect

The column $Type$ in table 7 indicates the type of attention distillation
used. $a,v,q,k$ represent the attention map, Value-Value relationship, Query-
Query relationship and Key-Key relationship respectively. From the
experimental results, it can be seen that the student network performance is
improved to some extent by adding attention map distillation. There was no
significant difference between distillation using only attention map and value
relations and distillation using attention graphs and all relations
simultaneously. The above results illustrate the effectiveness of attention
distillation.

To further investigate the effect of the parameter setting of $\alpha$ in Eq.
13 on the distillation effect, we set up a control experiment of attentional
distillation with different $\alpha$. The experiments use GRAPHSAGE, and the
distillation method uses the attention map and Value relations as additional
supervised signals. The experimental results are shown in Table 8 and Figure
4.

$\alpha$ weight | Distillation Imprv.
---|---
0.01 | +0.132
0.05 | +0.172
0.1 | +0.213
0.15 | +0.108
0.2 | +0.046
0.3 | +0.048
0.5 | +0.035
1 | -0.057
Table 8: :Relationship between $\alpha$ setting and distillation effect Figure
4: Relationship between $\alpha$ setting and distillation effect

It can be observed that the attention distillation effect increases with
increasing weight at the beginning and reaches its highest value at
$\alpha=0.1$. Subsequently, as the $\alpha=0.1$ increase, the distillation
effect begins to decline, even impairs the performance of student network.
This demonstrates that attention map distillation is better to be served as an
additional supervisory signal. By constraining the student network to emulate
the local topology of the teacher network during training, it allows the
student network to learn to better aggregate information from topology. When
this constraint is small, the student network can learn the appropriate node
representation according to the scale of its own features. When this
constraint is large, it has a negative impact on the learning of the student
network, where the student nodes are not free to explore the parameter space
and are forced to align to the distribution of the teacher network. However,
there is a scale inconsistency problem between the student model and the
teacher model. Forcing the student model to align to teacher model can hurt
the performance of the student model.

###  5.5 Comparison of GKEDM distillation and related work (RQ4)

In order to further demonstrate the effectiveness of GKEDM distillation, this
experiment compares the effectiveness of distillation using GKEDM with that of
related distillation methods. Other distillation methods specifically compared
are:

  * •

KD:It is the first method used to improve the performance of small models by
distilling the dark knowledge of a large teacher model. The core idea is to
use the output probability distribution of the teacher model as a soft label
to guide the learning process of the student model, thus improving the
generalization performance of the student model.

  * •

FITNET:FITNET uses an intermediate layer distillation approach. In the
distillation process, the intermediate layer output of the teacher model is
used as the target of the corresponding intermediate layer of the student
model, and the student model learns the knowledge of the teacher model by
minimizing the distance between its intermediate layer output and the
intermediate layer output of the teacher model. When the middle layer feature
size of the student network is different from that of the teacher network, the
student network output size and the teacher network output size can be made to
align by a linear mapping method.

  * •

LSP:It is the first method to apply knowledge distillation to graph neural
networks. Through the Local Structure Preservation (LSP) module, the teacher
network is able to guide the student network to learn the topology of graph
data, thus improving its performance.

The experimental results are shown in Table 9. In the KD method, the
hyperparameters of the student network training are set as follows: soft label
distillation weight is 0.8 and hard label weight is 0.2. The parameter
settings of the LSP are referred to the original paper. The weight of LSP loss
is 100, and the kernel method is set as RBF.

Model | Dataset |  Teacher Param.(M) | Teacher Acc. |  Student Param. (M) |  Student original Acc. | Distill Acc. | Type
---|---|---|---|---|---|---|---
GraphSAGE | PPI | 0.74 | 0.9213 | 0.15 | 0.8865 | 0.8527 | KD
GraphSAGE | PPI | 0.74 | 0.9213 | 0.15 | 0.8865 | 0.8823 | LSP
GraphSAGE | PPI | 0.74 | 0.9213 | 0.15 | 0.8865 | 0.8500 | FITNET
GraphSAGE | PPI | 0.74 | 0.9213 | 0.15 | 0.8865 | 0.9019 | a+v+q+k(ours)
GraphSAGE | FLICKR | 0.4 | 0.8183 | 0.2 | 0.6983 | 0.6977 | KD
GraphSAGE | FLICKR | 0.4 | 0.8183 | 0.2 | 0.6983 | 0.7049 | LSP
GraphSAGE | FLICKR | 0.4 | 0.8183 | 0.2 | 0.6983 | 0.6771 | FITNET
GraphSAGE | FLICKR | 0.4 | 0.8183 | 0.2 | 0.6983 | 0.7155 | a+v(ours)
MoNet | PPI | 1.04 | 0.9753 | 0.21 | 0.9023 | 0.9166 | KD
MoNet | PPI | 1.04 | 0.9753 | 0.21 | 0.9023 | 0.9297 | LSP
MoNet | PPI | 1.04 | 0.9753 | 0.21 | 0.9023 | 0.9282 | FITNET
MoNet | PPI | 1.04 | 0.9753 | 0.21 | 0.9023 | 0.9302 | a+v+q+k(ours)
Table 9: :Comparison of GKEDM distillation and related work

From the experimental results, it can be seen that attention distillation
using GKEDM achieved better results than the previous method on different GCNs
and different datasets, which illustrates the importance of attention for
improving the performance of GCNs. Compared with the unlearnable distance
measurement used in LSP, the distance measurement based on learnable
parameters of GKEDM is more flexible and reduces the constraints on the
student model, thus achieving better results.

##  6 Shortcomings and future work

Although GKEDM has achieved good results on various GCNs and datasets, the
following problems remain:

  * •

Since GKEDM is still a GCN implemented based on the message passing mechanism,
it also suffers from the problem of over-smoothing when the number of layers
is too large. It has been verified that the performance degradation still
occurs after stacking the attention layers of GKEDM. This problem reduces the
scope of GKEDM applications. The limitation of the number of parameters limits
the effectiveness of GKEDM on large datasets. The performance of GKEDM on
large-scale models and datasets has not been fully tested in this work.

  * •

Although GKEDM distillation is effective, the GCN has to be added with GKEDM
module and retrained. For large datasets, this process is very time and
resource consuming. In addition, the distillation process of GKEDM is
relatively complex and is not an end-to-end process.

Therefore, the future work could be:

  * •

Investigate the efficacy of GKEDM on large GCNs and large datasets, then
conceive a more universal framework.

  * •

Simplifying the distillation process, enables GKEDM to save computational time
and computational resources while maintaining the distillation effect.

##  7 Conclusion

In this work, we propose a knowledge enhancement and distillation module for
GCN, called GKEDM, which aims to improve the performance of GCNs on node
classification tasks and to lightweight them. The enhancement module of GKEDM
enhances the ability of GCNs to learn graph node representations by
introducing a convolutional layer based on an attention mechanism, thus
improving their performance on node classification tasks. Experiments show
that GKEDM can provide performance enhancement for different kinds of GCNs on
different datasets.

To make the model more lightweight through knowledge distillation, we also
propose a distillation module suitable for GKEDM. This module allows the
student network to mimic the neighborhood structure of the teacher network
through attention map distillation. The experimental results show that the
distillation module of GKEDM achieves the best results on several commonly
used node classification task datasets, proving the effectiveness of the
module. We believe that the method can provide useful guidance for broader
research on GCNs and can also serve as an effective solution to improve the
performance of node classification tasks. We believe that the method can
provide useful guidance for broader research on GCNs and can also serve as an
effective solution to improve the performance of GCNs on node classification
tasks. we hope this work can stimulate more research on GCNs and provide
support for broader applications in the future.

## References

  * [1] G. Li, M. Müller, B. Ghanem, and V. Koltun, “Training graph neural networks with 1000 layers,” Jun 2021.
  * [2] C. Ying, T. Cai, S. Luo, S. Zheng, G. Ke, D. He, Y. Shen, and T.-Y. Liu, “Do transformers really perform bad for graph representation,” Jun 2021.
  * [3] G. Hinton, O. Vinyals, and J. Dean, “Distilling the knowledge in a neural network,” Mar 2015.
  * [4] F. Tung and G. Mori, “Similarity-preserving knowledge distillation,” Jul 2019\.
  * [5] W. Park, D. Kim, Y. Lu, and M. Cho, “Relational knowledge distillation,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , Jan 2020. [Online]. Available: http://dx.doi.org/10.1109/cvpr.2019.00409
  * [6] Y. Yang, J. Qiu, M. Song, D. Tao, and X. Wang, “Distilling knowledge from graph convolutional networks,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , Aug 2020. [Online]. Available: http://dx.doi.org/10.1109/cvpr42600.2020.00710
  * [7] C. K. Joshi, F. Liu, X. Xun, J. Lin, and C. S. Foo, “On representation knowledge distillation for graph neural networks,” _IEEE Transactions on Neural Networks and Learning Systems_ , p. 1–12, Dec 2022. [Online]. Available: http://dx.doi.org/10.1109/tnnls.2022.3223018
  * [8] H. He, J. Wang, Z. Zhang, and F. Wu, “Compressing deep graph neural networks via adversarial knowledge distillation,” May 2022.
  * [9] T. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” Sep 2016.
  * [10] Z. Liu and J. Zhou, _Graph Attention Networks_ , Jun 2022, p. 39–41. [Online]. Available: http://dx.doi.org/10.1007/978-3-031-01587-8_7
  * [11] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” Jun 2017.
  * [12] M. Chen, Z. Wei, Z. Huang, B. Ding, and Y. Li, “Simple and deep graph convolutional networks,” Jul 2020.
  * [13] F. Monti, D. Boscaini, J. Masci, E. Rodola, J. Svoboda, and M. M. Bronstein, “Geometric deep learning on graphs and manifolds using mixture model cnns,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , Nov 2017. [Online]. Available: http://dx.doi.org/10.1109/cvpr.2017.576
  * [14] Y. Tian, D. Krishnan, and P. Isola, “Contrastive representation distillation,” Oct 2019.
  * [15] A. Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” Jul 2018.
  * [16] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” Jun 2017.
  * [17] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby, “An image is worth 16x16 words: Transformers for image recognition at scale,” Oct 2020.
  * [18] V. Dwivedi, C. Joshi, T. Laurent, Y. Bengio, and X. Bresson, “Benchmarking graph neural networks,” Mar 2020.
  * [19] W. Wang, F. Wei, L. Dong, H. Bao, N. Yang, and M. Zhou, “Minilm: Deep self-attention distillation for task-agnostic compression of pre-trained transformers,” Feb 2020.
  * [20] M. Zitnik and J. Leskovec, “Predicting multicellular function through multi-layer tissue networks.” _Bioinformatics_ , vol. 33, no. 14, p. i190–i198, Apr 2017. [Online]. Available: http://dx.doi.org/10.1093/bioinformatics/btx252
  * [21] H. Zeng, H. Zhou, A. Srivastava, R. Kannan, and V. Prasanna, “Graphsaint: Graph sampling based inductive learning method,” Jul 2019.
  * [22] O. Shchur, M. Mumme, A. Bojchevski, and S. Günnemann, “Pitfalls of graph neural network evaluation,” _Relational Representation Learning Workshop, NeurIPS 2018_ , 2018.


e1e-mail<EMAIL_ADDRESS>11institutetext: Department of
Mechanical Science and Bioengineering, Graduate School of Engineering Science,
Osaka University, 1–3 Machikaneyama, Toyonaka, Osaka 560–8531, Japan

We numerically investigate the mechanical and geometrical properties of dense
wet granular particles with irreversible attractive interaction. The shear
modulus exhibits two inflection points as the packing fraction increases, and
the bulk modulus shows a non-monotonic behavior. The coordination number also
exhibits two inflection points. The peak position in the pair correlation
function shifts to a lower value due to attractive interaction. The Voronoi
tessellation of the particle configuration reveals that the probability
density function for the volume of the Voronoi cell broadens as the packing
fraction approaches the jamming point.

# Mechanical and geometrical properties of jammed wet granular materials

Kiwamu Yoshiie1,addr1 and Michio Otsukiaddr1

(Received: date / Accepted: date)

##  1 Introduction

Amorphous materials, such as granular materials, emulsions, foams, and
colloidal suspensions, behave like solids with rigidity when the packing
fraction $\phi$ exceeds a critical volume fraction $\phi_{c}$, while they
exhibit liquid-like behaviors for $\phi<\phi_{c}$. This transition, known as
the jamming transition, has been extensively studied for years liu1998jamming
; liu2010jamming ; van2009jamming ; behringer2018physics .

The mechanical properties of the materials exhibit critical behaviors near the
critical fraction $\phi_{c}$. The pressure $P$, the shear modulus $G$, and the
bulk modulus $B$ for repulsive frictionless particles exhibit power law
scalings as a function of $\phi-\phi_{c}$ ohern2002random ; ohern2003jamming .
Critical scaling laws for the rheological properties have been observed for
systems under steady shear olsson2007critical ; hatano2007criticality ;
hatano2008scaling ; tighe2010model ; otsuki2009critical . Recently, the
nonlinear elasticity of jammed amorphous materials coulais2014shear ;
otsuki2014avalanche ; boschan2016beyond ; otsuki2022softening , the effect of
friction somfai2007critical ; silbert2010jamming ; otsuki2017discontinuous ;
otsuki2021shear , and the frequency dependence of the complex shear modulus
tighe2011relaxations ; dagois2017softening  have been studied.

The geometrical properties also change drastically in the vicinity of
$\phi_{c}$. When the packing fraction exceeds $\phi_{c}$, the coordination
number $Z$ of frictionless particles changes from zero to the isostatic value
$Z_{\mathrm{iso}}$. The excess coordination number $Z-Z_{\mathrm{iso}}$
exhibits a power law scaling as a function of $\phi-\phi_{c}$ ohern2002random
; ohern2003jamming . Moreover, in the three-dimensional systems consisting of
frictionless mono-dispersed spheres with diameter $d$, the pair correlation
function $g(r)$ diverges as $r\to d$ at $\phi_{c}$, which indicates that a lot
of particles are on the verge of making contact silbert2006structural .

Earlier studies on the jamming transition have focused on purely repulsive
particles. However, cohesion, such as the capillary force in wet granular
materials, is nonnegligible in various realistic situations. It is well known
that the irreversible capillary force affects the dynamics of granular
materials herminghaus2005dynamics ; strauch2012wet ; herminghaus2013wet ;
mitarai2006wet , which might result from the change in their mechanical and
geometrical properties. Some researchers have studied the rheological
properties of cohesive particles in systems at constant volume
chaudhuri2012inhomogeneous ; gu2014rheology ; irani2014impact ;
irani2016athermal ; irani2019discontinuous  or under constant pressure
rognon2008dense ; khamseh2015flow ; yamaguchi2018rheology ; badetti2018shear ;
mandal2021rheology ; vo2020additive ; vo2020evolution ; vo2020role ;
macaulay2021viscosity . The yield stress or the apparent friction coefficient
increases due to the attractive interaction between particles. The attractive
force also causes a gel-like contact network for low $\phi$ head2007well ;
zheng2016shear , shear bands irani2014impact ; irani2016athermal ;
irani2019discontinuous ; singh2014effect , and clusters of particles
yamaguchi2018rheology ; vo2020evolution ; macaulay2019shear ; lois2008jamming
. For two-dimensional particles with a simple reversible attractive force, it
is reported that the transition point $\phi_{c}$ for rigidity decreases as the
strength of cohesion increases koeze2018sticky , and the critical scaling laws
in $G$ and $B$ for repulsive particles are broken Koeze2020Elasticity .
However, the behavior of three-dimensional wet granular materials with the
irreversible capillary force near $\phi_{c}$ remains unclear.

In this study, we numerically investigate the mechanical and geometrical
properties of three-dimensional frictionless wet granular materials. In Sect.
2, we explain our model and setup. Section 3 presents the numerical results on
the mechanical properties. In Sect. 4, we deal with the geometrical
properties. Section 4 consists of three parts. The $\phi$-dependence of the
coordination number is shown in Sect. 4.1. We demonstrate the change in the
pair correlation function for $\Gamma_{s}>0$ in Sect. 4.2. In Sect. 4.3, we
analyze the geometry of the particle configuration using the Voronoi
tessellation. Section 5 presents discussion and conclusions of our results.
The dependence of the critical fraction $\phi_{c}$ on the strength of cohesion
is shown in A. In B, we investigate the effect of the capillary force on
scaling laws for $G$.

##  2 Setup

We consider three-dimensional wet granular materials consisting of $N$
frictionless identical particles in a cubic box of linear size $L$. We apply
oscillatory shear under Lees-Edwards boundary conditions using the SLLOD
method to measure the mechanical properties evans2008non . The equation of
motion is given by

$\displaystyle\dfrac{\differential\bm{r}_{i}}{\differential t}$
$\displaystyle=\dfrac{\bm{p}_{i}}{m}+\dot{\gamma}(t)y_{i}\bm{e}_{x},$ (1)
$\displaystyle\dfrac{\differential\bm{p}_{i}}{\differential t}$
$\displaystyle=\sum_{j\neq
i}F_{ij}\bm{n}_{ij}-\dot{\gamma}(t)p_{i,y}\bm{e}_{x}$ (2)

with the mass $m$, the unit vector along the x-direction $\bm{e}_{x}$, and the
shear rate $\dot{\gamma}(t)$. Here, $\bm{r}_{i}=(x_{i},y_{i},z_{i})$ and
$\bm{p}_{i}=m_{i}(\frac{d}{dt}{\bm{r}}_{i}-\dot{\gamma}(t)y_{i})\bm{e}_{x}$
are the position and peculiar momentum of particle $i$, respectively. The
force between particles $i$ and $j$ is denoted by $F_{ij}$. The force $F_{ij}$
consists of the dissipative force $F_{ij}^{(\mathrm{d})}$ and the static force
$F_{ij}^{(\mathrm{s})}$ as

$\displaystyle F_{ij}=F_{ij}^{(\mathrm{d})}+F_{ij}^{(\mathrm{s})}.$ (3)

The normal unit vector $\bm{n}_{ij}$ is given by
$\bm{n}_{ij}=\bm{r}_{ij}/r_{ij}$ with $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ and
$r_{ij}=|\bm{r}_{ij}|$.

The dissipative force is given by

$\displaystyle F_{ij}^{(\mathrm{d})}$
$\displaystyle=-\xi_{\mathrm{n}}v^{\rm(n)}_{ij}\Theta(\Delta_{ij})$ (4)

with the normal velocity
$v_{ij}^{\rm(n)}=(\frac{d}{dt}\bm{r}_{i}-\frac{d}{dt}\bm{r}_{j})\cdot\bm{n}_{ij}$
and the viscous constant $\xi_{\mathrm{n}}$, the diameter $d$, and the overlap
length $\Delta_{ij}=d-r_{ij}$. The contact between particles is formed at
$\Delta_{ij}=0$ ($r_{ij}=d$). Here, $\Theta(x)$ is the Heaviside step function
satisfying $\Theta(x)=1$ for $x\geq 0$ and $\Theta(x)=0$ otherwise. The static
force $F_{ij}^{(\mathrm{s})}$ consists of the elastic contact force
$F_{ij}^{(\mathrm{e})}$ and the attractive capillary force
$F_{ij}^{(\mathrm{cap})}$ as

$\displaystyle F_{ij}^{(\mathrm{s})}$
$\displaystyle=F_{ij}^{(\mathrm{e})}+F_{ij}^{(\mathrm{cap})}.$ (5)

The elastic contact force is given by

$\displaystyle F_{ij}^{(\mathrm{e})}$
$\displaystyle=k_{\mathrm{n}}\Delta_{ij}\Theta(\Delta_{ij})$ (6)

with the spring constant $k_{\mathrm{n}}$.

Figure 1: Static force $F_{ij}^{(\mathrm{s})}$ against $\Delta_{ij}$. The red
line represents the behavior when the particles are approaching. The blue line
represents the behavior when the particles are separating.

The capillary force $F_{ij}^{(\mathrm{cap})}$ is modeled as roy2016micro ;
roy2017general ; roy2018liquid ; shi2020steady

$\displaystyle
F_{ij}^{(\mathrm{cap})}=\begin{cases}-2\pi\Gamma_{s}d\cos\theta&\mbox{if
}\Delta_{ij}\geq 0,\\\\[4.0pt]
f^{\rm(a)}(\hat{s}_{ij})\Theta\left(\Delta_{ij}+d_{\mathrm{c}}\right)&\mbox{if
}\Delta_{ij}<0,\mbox{separation},\\\\[4.0pt] 0&\mbox{if
}\Delta_{ij}<0,\mbox{approaching}\\\\[4.0pt] \end{cases}$ (7)

with

$\displaystyle
f^{\rm(a)}(\hat{s}_{ij})=\dfrac{-2\pi\Gamma_{s}d\cos\theta}{1+1.05\hat{s}_{ij}+2.5\hat{s}_{ij}^{2}}.$
(8)

Here, $\theta$ is the contact angle, $d_{c}$ is the rapture length of the
capillary bridge, $\hat{s}_{ij}=s_{ij}\sqrt{d/V_{b}}$ is the normalized
separation distance, with the separation distance $s_{ij}=-\Delta_{ij}$. The
surface tension is denoted by $\Gamma_{s}$, which characterizes the strength
of cohesion. The volume of the liquid bridge is denoted by $V_{b}$, which
characterizes the rapture length $d_{c}$ as $d_{c}=(1+\theta/2)V_{b}^{1/3}$
lian1993theoretical ; willett2000capillary . When the particles approach
before the contact, $F_{ij}^{(\mathrm{cap})}$ is zero. After the particles
contact at $\Delta_{ij}=0$, the capillary bridge is formed, i.e., the
attractive force becomes active. When the particles are separating, the force
does not follow the same path; the attractive force is active until the
capillary bridge disappears at $\Delta_{ij}=-d_{c}$.

The static force $F^{(\mathrm{s})}$ given by Eqs. (5)-(8) is irreversible due
to the capillary force, as shown in Fig. 1. When the particles are approaching
before contact, $F^{(\mathrm{s})}$ is zero and follows the red line. After the
capillary bridge is formed at $\Delta_{ij}=0$, the capillary force
$F^{(\mathrm{cap})}$ is active, and $F^{(\mathrm{s})}$ follows the blue line
until the bridge is broken at $\Delta_{ij}=-d_{c}$. As shown in Fig. 1, the
static force becomes zero at $\Delta_{ij}=\delta_{0}$ ($r_{ij}=d-\delta_{0}$)
with $\delta_{0}=2\pi\Gamma_{s}d\cos\theta/k_{\mathrm{n}}>0$, where positive
$F_{ij}^{(\mathrm{e})}$ and negative $F_{ij}^{(\mathrm{cap})}$ cancel each
other out.

The particles are randomly placed with an initial packing fraction
$\phi_{\mathrm{ini}}=0.45$ without any overlap. The system is gradually
compressed until the packing fraction reaches a given value $\phi$. In each
compression step, we increase the packing fraction by $\Delta\phi=0.000025$
with the affine transformation of the particle configuration and the system
size. The particles are relaxed to a mechanical equilibrium state with the
kinetic temperature $T<T_{\mathrm{th}}$ following Eqs. (1) and (2) with
$\dot{\gamma}(t)=0$. Here, the kinetic temperature is given by $T=\sum
m|\bm{v}_{i}|^{2}/(2N)$.

After the compression, we apply the oscillatory shear strain as

$\displaystyle\gamma(t)=\gamma_{0}\sin\omega t$ (9)

for $N_{\mathrm{cyc}}$ cycles. Here, $\gamma_{0}$ and $\omega$ are the strain
amplitude and the angular frequency, respectively, which are set small enough.
In the last cycle, we measure the shear (storage) modulus $G$ as doi1988theory

$\displaystyle G(\phi)=\frac{\omega}{\pi}\int_{0}^{2\pi/\omega}\differential
t\ \sigma_{xy}(\phi,t)\sin\omega t/\gamma_{0}$ (10)

with the shear stress for a given $\phi$

$\displaystyle\sigma_{xy}(\phi,t)=-\dfrac{1}{2L^{3}}\sum_{i}\sum_{i<j}(r_{ij,x}F_{ij,y}+r_{ij,y}F_{ij,x}).$
(11)

We also measure the pressure as

$\displaystyle
P(\phi)=\dfrac{1}{3L^{3}}\sum_{i}\sum_{i<j}\bm{r}_{ij}\cdot\bm{F}_{ij}$ (12)

after the last cycle, and calculate the bulk modulus as

$\displaystyle B(\phi)=\phi\frac{\differential P}{\differential\phi}.$ (13)

Here, we have ignored the kinetic parts of the shear stress and the pressure
because the contact stress is dominant in our dense system da2005rheophysics .

We use $N=3000$, $\gamma_{0}=1.0\times 10^{-5}$, $\omega=1.0\times
10^{-4}\sqrt{k_{\mathrm{n}}/m}$, $N_{\mathrm{cyc}}=100$, and
$T_{\mathrm{th}}=10^{-8}k_{\mathrm{n}}d_{0}^{2}$. We choose $d_{c}=5.0\times
10^{-4}d_{0}$, $V_{b}=7.5\times 10^{-11}d_{0}^{3}$, $\theta=\pi/9$, and
$\Gamma_{s}/k_{\mathrm{n}}=0,~{}3.0\times 10^{-3},~{}1.5\times
10^{-2},~{}3.0\times 10^{-2}$ for the parameters of the capillary force
following ref. roy2017general ; roy2018liquid ; shi2020steady . Here,
$\Gamma_{s}=0$ corresponds to dry particles. We adopt the Adams-Morton and
Adams-Bashforth methods with a time step $\Delta
t=0.005\sqrt{m/k_{\mathrm{n}}}$ for the time evolution of $\bm{r}_{i}$ and
$\bm{p}_{i}$, respectively. We have numerically confirmed that $N$ and
$N_{\mathrm{cyc}}$ are large enough, and $\gamma_{0}$, $\omega$,
$T_{\mathrm{th}}$, and $\Delta t$ are small enough not to influence our
results.

##  3 Mechanical properties

In Fig. 2, we plot the shear modulus $G$ against the volume fraction $\phi$
for various $\Gamma_{s}$. For each $\Gamma_{s}$, $G$ becomes non-zero as the
packing fraction $\phi$ exceeds a critical fraction $\phi_{c}$. As
$\Gamma_{s}$ increases, the critical fraction $\phi_{c}$ decreases, as shown
in A. The shear modulus $G$ increases with $\phi$. For wet particles with
$\Gamma_{s}>0$, there are two inflection points (open symbols) where the
curvature of $G(\phi)$ changes sign with
$\frac{\differential[2]}{\differential\phi^{2}}G(\phi)=0$. The position of the
inflection point with higher $\phi$ is almost independent of $\Gamma_{s}$
($\phi\simeq 0.63$). The inflection point does not exist for dry particles
with $\Gamma_{s}=0$. The inflection points indicate that the simple power law
scaling for repulsive particles ohern2002random ; ohern2003jamming  is not
satisfied for wet granular materials. Similar behaviors are reported for two-
dimensional particles with a simple reversible attractive interaction
Koeze2020Elasticity , but only one inflection point exists in the system.

Figure 2: Shear modulus $G$ against $\phi$ for various $\Gamma_{s}$ with
$\phi>\phi_{c}$. Open symbols represent inflection points. Figure 3: Pressure
$P$ against $\phi$ for various $\Gamma_{s}$ with $\phi>\phi_{c}$. The inset
shows $P$ in the vicinity of $\phi_{c}$. The dashed line represents $P=0$.

Figure 3 displays the pressure $P$ against the volume fraction $\phi$ for
various $\Gamma_{s}$ with $\phi>\phi_{c}$. The inset of Fig. 3 shows $P$ near
$\phi_{c}$. For each $\Gamma_{s}$, $P$ increases with $\phi$. The pressure is
positive and almost $0$ even in the vicinity of $\phi_{c}$, as shown in the
inset of Fig. 3. For high $\phi>0.65$, $P$ seems independent of $\Gamma_{s}$.
See B for the relation between $G$ and $P$.

Figure 4: Bulk modulus $B$ against $\phi$ for various $\Gamma_{s}$ with
$\phi>\phi_{c}$.

In Fig. 4, we demonstrate the bulk modulus $B$ against the packing fraction
$\phi$ for various $\Gamma_{s}$ with $\phi>\phi_{c}$. The bulk modulus $B$ is
not a monotonic function of $\phi$ for $\Gamma_{s}>0$, while $B$ for
$\Gamma_{s}=0$ does not exhibit such non-monotonic behavior. As the packing
fraction $\phi$ decreases, $B$ rapidly increases near $\phi_{c}$, which is not
shown in the previous study Koeze2020Elasticity . The non-monotonic behavior
in $B$ of wet granular materials indicates that $B$ does not obey the power
law scaling for dry repulsive particles ohern2002random ; ohern2003jamming .

##  4 Geometrical properties

In this section, we analyze the geometrical properties of wet granular
materials. In Sect. 4.1, we show the $\phi$-dependence of the coordination
number. Section 4.2 demonstrates the change in the pair correlation function
for $\Gamma_{s}>0$. In Sect. 4.3, we discuss the probability density function
for the volume of the Voronoi cell obtained by the Voronoi tessellation.

###  4.1 Coordination number

We plot the excess coordination number $Z-Z_{\mathrm{iso}}$ against the volume
fraction $\phi$ for various $\Gamma_{s}$ with $\phi>\phi_{c}$ in Fig. 5. Here,
we calculate the coordination number $Z$ after the final oscillatory shear as

$\displaystyle Z=2N_{\mathrm{con}}/N,$ (14)

where $N_{\mathrm{con}}$ is the total number of contacts with $F_{ij}\neq 0$.
For three-dimensional frictionless particles, the isostatic value
$Z_{\mathrm{iso}}$ equals $6$ van2009jamming . As $\phi$ increases, so does
$Z-Z_{\mathrm{iso}}$, increasing from $0$. There are two inflection points
(open symbols) in $Z-Z_{\mathrm{iso}}$ for each $\Gamma_{s}$. Their positions
are almost the same as those for $G$ in Fig. 2. The existence of two
inflection points is natural if the relation $G\propto Z-Z_{\mathrm{iso}}$ for
repulsive particles ohern2002random  holds for wet granular particles, which
is discussed in B.

Figure 5: Excess coordination number $Z-Z_{\mathrm{iso}}$ against $\phi$ for
various $\Gamma_{s}$ with $\phi>\phi_{c}$. Open symbols represent the
inflection points.

###  4.2 Pair correlation function

Figure 6: Pair correlation function $g(r)$ against $r$ for various $\phi$ with
$\Gamma_{s}/k_{\mathrm{n}}=0$ (a) and $\Gamma_{s}/k_{\mathrm{n}}=3.0\times
10^{-3}$ (b). Dashed and dash-dotted lines represent $r=d$ and
$r=d-\delta_{0}$, respectively. The blue shaded region corresponds to
$r<d-\delta_{0}$. The red shaded region corresponds to $d-\delta_{0}\leq r\leq
d+d_{c}$.  Figure 7: Plot of $\Psi_{\mathrm{+}}$ and $\Psi_{\mathrm{-}}$
against $\phi$ for various $\Gamma_{s}$ with $\phi>\phi_{c}$. Circles ,
squares, and triangles represent $\Gamma_{s}=3.0\times 10^{-3}$,
$\Gamma_{s}=1.5\times 10^{-3}$, and $\Gamma_{s}=3.0\times 10^{-4}$
respectively. Open (filled) symbol corresponds to $\Psi_{\mathrm{+}}$
($\Psi_{\mathrm{-}}$).

In Fig. 6, we demonstrate the pair correlation function $g(r)$ against $r$ in
the static state after the final oscillatory shear for various
$\phi>\phi_{c}$. The pair correlation function $g(r)$ is given by

$\displaystyle g(r)=\dfrac{L^{3}}{N^{2}}\left\langle\sum_{i}\sum_{j\neq
i}\delta^{3}(r-r_{ij})\right\rangle.$ (15)

The blue area in Fig. 6 ($r<d-\delta_{0}$) corresponds to the region with the
static force $F_{ij}^{(\mathrm{s})}>0$, while the red area
($d-\delta_{0}<r<d+d_{c}$) represents the region with
$F_{ij}^{(\mathrm{s})}<0$. The first peak of $g(r)$ rapidly increases as
$\phi\to\phi_{c}$. For $\Gamma_{s}/k_{\mathrm{n}}=0$ (Fig. 6 (a)), the peak
position approaches $r=d$ as $\phi\to\phi_{c}$ silbert2006structural . For
$\Gamma_{s}/k_{\mathrm{n}}=3.0\times 10^{-3}$ (Fig. 6 (b)), the peak position
decreases to $r=d-\delta_{0}$, at which $F_{ij}^{(\rm s)}=0$. We should note
that the gel-like structure characterized by power law decay in the structure
factor zheng2016shear  is not observed in our simulation with $\phi>\phi_{c}$.

In Fig. 6 (b), $g(r)$ has a large value even in the red region with
$d-\delta_{0}<r<d+d_{c}$ near $\phi_{c}$, which indicates that many contacts
have a negative static force $F_{ij}^{(\mathrm{s})}$. Here, we introduce
$N_{+}$ and $N_{-}$ as the number of contacts with $F_{ij}^{(\mathrm{s})}>0$
and $F_{ij}^{(\mathrm{s})}<0$, respectively, and plot the ratios
$\Psi_{+}=N_{+}/N_{\mathrm{con}}$ and $\Psi_{-}=N_{-}/N_{\mathrm{con}}$ with
the total number of contacts $N_{\mathrm{con}}=N_{+}+N_{-}$ against $\phi$ for
various $\Gamma_{s}>0$ in Fig. 7. Note that
$N_{+}=\int_{0}^{d-\delta_{0}}g(r)\differential r$ and
$N_{-}=\int^{d+d_{c}}_{d-\delta_{0}}g(r)\differential r$ if all pairs with
$r_{ij}<d-\delta_{0}$ are contacting. For high $\phi$, $F_{ij}^{(\mathrm{s})}$
is positive for almost all contacts ($\Psi_{+}=1$ and $\Psi_{-}=0$), which is
consistent with high $P$ for large $\phi$ in Fig. 3. As $\phi$ decreases to
$\phi_{c}$, $\Psi_{+}$ and $\Psi_{-}$ approach $0.5$. This indicates that the
pressure $P$ of wet granular materials decreases as $\phi\to\phi_{c}$ because
the number of contacts with $F_{ij}^{(\mathrm{s})}<0$ increases. This behavior
is different from that of dry particles with purely repulsive interaction,
where the positive contact force decreases to $0$, keeping
$N_{+}=N_{\mathrm{con}}\simeq Z_{\rm iso}N/2$ ($\Psi_{+}=1$).

###  4.3 Voronoi tessellation

The structure of disordered particles has been studied using the Voronoi
tessellation bernal1964bakerian ; finney1970random ; finney1970random2 ;
oger1996voronoi ; jullien1996computer ; yang2002voronoi ; aste2007invariant ;
xu2007analysis . A Voronoi cell associated with each particle contains an
ensemble of points closer to a given sphere center than any other. In ref.
xu2007analysis , it is reported that the probability density function for the
volume of the Voronoi cell changes due to the capillary force, but the
$\phi$-dependence is not investigated.

Figure 8: Probability density function $p(V^{*})$ against $V^{*}$ with
$\Gamma_{s}/k_{\mathrm{n}}=0$ (a) and $\Gamma_{s}/k_{\mathrm{n}}=3.0\times
10^{-3}$ (b) for various $\phi>\phi_{c}$.

We obtain the Voronoi cell from the particle configuration after the final
oscillatory shear using the VORO++ code library rycroft2009voro++ . We define
the volume of the Voronoi cell as $V$ and plot the probability density
functions $p(V^{*})$ of the normalized volume $V^{*}=V/\bar{V}$ with the
average of the volume $\bar{V}=L^{3}/N$ in Fig. 8 for different
$\phi>\phi_{c}$. Note that $\bar{V}$ is related to the packing fraction $\phi$
as $\bar{V}=\pi d^{3}/(6\phi)$. For dry particles
($\Gamma_{s}/k_{\mathrm{n}}=0$), the probability density function is almost
independent of $\phi$ (Fig. 8(a)). This independence might indicate that the
volume of the Voronoi cell is affinely deformed by changing $\phi$. For wet
particles ($\Gamma_{s}/k_{\mathrm{n}}=3.0\times 10^{-3}$), the width of the
probability density functions increases, and the peak position of $p(V^{*})$
is shifted to lower $V^{*}$ as $\phi$ decreases (Fig. 8(b)), which is
consistent with an experiment of wet particles xu2007analysis .

Figure 9: Variance of $V^{*}$ against $\phi$ for various $\Gamma_{s}$ with
$\phi>\phi_{c}$.

In Fig. 9, the variance of the normalized volume $V^{*}$ is plotted against
$\phi$ for different $\Gamma_{s}$ with $\phi>\phi_{c}$. The variances for
different $\Gamma_{s}$ collapse onto a master curve. However, the range of
$\phi$ for each $\Gamma_{s}$ depends on $\phi_{c}$ and broadens as
$\Gamma_{s}$ increases. For $\Gamma_{s}/k_{\mathrm{n}}=0$, the variance is
almost independent of $\phi$, which is consistent with the probability density
function of $V^{*}$ in Fig. 8 (a). However, the variance for $\Gamma_{s}>0$
increases as $\phi$ approaches $\phi_{c}$. This corresponds to the increase in
the width of $p(V^{*})$ in Fig. 8 (b). The packing fraction where the variance
increases is close to the inflection point for $G$ ($\phi\simeq 0.63$) in Fig.
2.

##  5 Conclusion and discussion

We numerically studied the mechanical and geometrical properties of wet
granular materials with $\phi>\phi_{c}$. For $\Gamma_{s}>0$, the shear modulus
$G(\phi)$ has two inflection points, and the bulk modulus $B(\phi)$ exhibits a
non-monotonic behavior. These mechanical properties are qualitatively
different from those for dry particles with $\Gamma_{s}=0$, where $G(\phi)$
and $B(\phi)$ obey simple power law scalings near $\phi_{c}$ ohern2002random ;
ohern2003jamming . The excess coordination number $Z(\phi)-Z_{\mathrm{iso}}$
also has two inflection points. The peak position in the pair correlation
function $g(r)$ becomes lower than the diameter $d$ due to the attractive
capillary force. The probability density function for the volume of the
Voronoi cell broadens as the packing fraction approaches $\phi_{c}$. These
results indicate that the geometrical properties change with the mechanical
properties due to the capillary force.

The breakdown of simple power-law behaviors in $G(\phi)$ and $B(\phi)$ has
been reported in ref. Koeze2020Elasticity  for two-dimensional particles with
a simple reversible attractive force. However, it is also reported that $G$
and $B$ satisfy critical scaling laws, including the strength of attraction.
The attractive force in ref. Koeze2020Elasticity  differs from the
irreversible capillary force, and the critical scaling laws do not apply to
the mechanical properties shown in this paper because of the two inflection
points for $G(\phi)$ and the non-monotonic behavior in $B(\phi)$, which are
not observed in ref. Koeze2020Elasticity . The critical scaling laws near
$\phi_{c}$ in wet granular materials will be the subject of future study.

In this study, we have neglected contact friction between particles to focus
on the effect of the attractive capillary force. Recent studies have reported
that the contact friction affects the critical behaviors near $\phi_{c}$ for
dry repulsive particles somfai2007critical ; silbert2010jamming ;
otsuki2017discontinuous ; otsuki2021shear . However, it is unclear whether the
friction force changes the mechanical properties of attractive wet particles
shown in this study. Further work is necessary to resolve this issue.

## Acknowledgment

K. Y. and M.O. thank S. Takada, T. Nakamura, and H. Mizuno for helpful
discussions. Numerical computation in this work was conducted at the Yukawa
Institute Computer Facility. We would like to thank Editage (www.editage.com)
for English language editing. K.Y. is partially supported by Leave a Nest Co.,
Ltd., Hosokawa Powder Technology Foundation (Grant No. HPTF20506), and the
Grant-in-Aid for Japan Society for Promotion of Science JSPS Research Fellow
(Grant No. 21J13720). M.O. is partially supported by Scientific Grant-in-Aid
of Japan Society for the Promotion of Science, KAKENHI (Grants No. 19K03670
and No. 21H01006).

## Author contribution statement

K.Y. carried out the numerical simulations. K.Y. and M.O. interpreted the
results and wrote the manuscript.

##  Appendix A Critical fraction

This appendix shows the $\Gamma_{s}$-dependence of the critical fraction
$\phi_{c}$. Here, we define the critical fraction $\phi_{c}$ as the packing
fraction where $G$ exceeds a threshold $G_{\mathrm{th}}$ with
$G_{\mathrm{th}}/k_{\mathrm{n}}=1.0\times 10^{-4}$. We have checked that
$\phi_{c}$ does not change if we use a smaller
$G_{\mathrm{th}}/k_{\mathrm{n}}=5.0\times 10^{-5}$. In Fig. 10, we plot the
critical fraction $\phi_{c}$ against $\Gamma_{s}$. The critical fraction
$\phi_{c}$ decreases as $\Gamma_{s}$ increases, which is consistent with the
results of ref. Koeze2020Elasticity .

Figure 10: Critical fraction $\phi_{c}$ against $\Gamma_{s}$.

##  Appendix B Critical scaling of $G$

Figure 11: Shear modulus $G$ against $P$ for various $\Gamma_{s}$. The dashed
line represents $G\propto P^{1/2}$.

In this appendix, we numerically investigate the effect of the capillary force
on the scaling laws for $G$ of dry repulsive particles. For frictionless
particles with the linear elastic interaction $F_{ij}^{\mathrm{(e)}}$ given by
Eq. (6), $G$ satisfies $G\propto P^{\alpha}$ with $\alpha\simeq 1/2$ and
$G\propto(Z-Z_{\mathrm{iso}})$ ohern2002random ; ohern2003jamming .

Figure 11 displays $G$ against $P$ for various $\Gamma_{s}$ obtained from the
data in Figs. 2 and 3. For dry particles with $\Gamma_{s}=0$, $G$ almost
satisfies the scaling law $G\propto P^{\alpha}$ with $\alpha=1/2$. For wet
particles with $\Gamma_{s}>0$, there is a region where $G$ behaves as a power
law function with an exponent lower than $1/2$ for
$P/(k_{\mathrm{n}}d^{-1})>10^{4}$. However, the power-law behavior seems to
break down as $P$ decreases, and $G$ seems to become zero at a finite $P$.

In Fig. 12, we show the shear modulus $G$ against the excess coordination
number $Z-Z_{\mathrm{iso}}$ for different $\Gamma_{s}$ obtained from the data
in Figs. 2 and 5. For dry particles with $\Gamma_{s}=0$,
$G\propto(Z-Z_{\mathrm{iso}})$ is satisfied. For $\Gamma_{s}>0$, there is a
region where $G$ is proportional to $Z-Z_{\mathrm{iso}}$, but the scaling
relation is broken for smaller $Z-Z_{\mathrm{iso}}$ near $\phi_{c}$.

Figure 12: Shear modulus $G$ against $Z-Z_{\mathrm{iso}}$ for various
$\Gamma_{s}$. The dashed line represents $G\propto Z-Z_{\mathrm{iso}}$.

## References

  * (1) A.J. Liu, S.R. Nagel, Nature 396(6706), 21 (1998)
  * (2) A.J. Liu, S.R. Nagel, Annu. Rev. Condens. Matter Phys 1(1), 347 (2010)
  * (3) M. van Hecke, J. Phys. Condens. Matter 22(3), 033101 (2009)
  * (4) R.P. Behringer, B. Chakraborty, Rep. Prog. Phys. 82(1), 012601 (2019)
  * (5) C.S. O’Hern, S.A. Langer, A.J. Liu, S.R. Nagel, Phys. Rev. Lett. 88(7), 075507 (2002)
  * (6) C.S. O’Hern, L.E. Silbert, A.J. Liu, S.R. Nagel, Phys. Rev. E 68(1), 011306 (2003)
  * (7) P. Olsson, S. Teitel, Phys. Rev. Lett. 99(17), 178001 (2007)
  * (8) T. Hatano, M. Otsuki, S.i. Sasa, J. Phys. Soc. Jpn. 76(2), 023001 (2007)
  * (9) T. Hatano, J. Phys. Soc. Jpn. 77(12), 123002 (2008)
  * (10) B.P. Tighe, E. Woldhuis, J.J. Remmers, W. van Saarloos, M. van Hecke, Phys. Rev. Lett. 105(8), 088303 (2010)
  * (11) M. Otsuki, H. Hayakawa, Phys. Rev. E 80(1), 011308 (2009)
  * (12) C. Coulais, A. Seguin, O. Dauchot, Phys. Rev. Lett. 113(19), 198001 (2014)
  * (13) M. Otsuki, H. Hayakawa, Phys. Rev. E 90(4), 042202 (2014)
  * (14) J. Boschan, D. Vågberg, E. Somfai, B.P. Tighe, Soft Matter 12(24), 5450 (2016)
  * (15) M. Otsuki, H. Hayakawa, Phys. Rev. Lett. 128(20), 208002 (2022)
  * (16) E. Somfai, M. van Hecke, W.G. Ellenbroek, K. Shundyak, W. van Saarloos, Phys. Rev. E 75(2), 020301 (2007)
  * (17) L.E. Silbert, Soft Matter 6(13), 2918 (2010)
  * (18) M. Otsuki, H. Hayakawa, Phys. Rev. E 95(6), 062902 (2017)
  * (19) M. Otsuki, H. Hayakawa, Eur. Phys. J. E 44(5), 1 (2021)
  * (20) B.P. Tighe, Phys. Rev. Lett. 107(15), 158303 (2011)
  * (21) S. Dagois-Bohy, E. Somfai, B.P. Tighe, M. van Hecke, Soft matter 13(47), 9036 (2017)
  * (22) L.E. Silbert, A.J. Liu, S.R. Nagel, Phys. Rev. E 73(4), 041304 (2006)
  * (23) S. Herminghaus, Adv. Phys. 54(3), 221 (2005)
  * (24) S. Strauch, S. Herminghaus, Soft Matter 8(32), 8271 (2012)
  * (25) S. Herminghaus, _Wet granular matter: a truly complex fluid_ , vol. 6 (World Scientific, 2013)
  * (26) N. Mitarai, F. Nori, Adv. Phys. 55(1-2), 1 (2006)
  * (27) P. Chaudhuri, L. Berthier, L. Bocquet, Phys. Rev. E 85(2), 021503 (2012)
  * (28) Y. Gu, S. Chialvo, S. Sundaresan, Phys. Rev. E 90(3), 032206 (2014)
  * (29) E. Irani, P. Chaudhuri, C. Heussinger, Phys. Rev. Lett. 112(18), 188303 (2014)
  * (30) E. Irani, P. Chaudhuri, C. Heussinger, Phys. Rev. E 94(5), 052608 (2016)
  * (31) E. Irani, P. Chaudhuri, C. Heussinger, Phys. Rev. Fluids 4(7), 074307 (2019)
  * (32) P.G. Rognon, J.N. Roux, M. Naaim, F. Chevoir, J. Fluid Mech. 596, 21 (2008)
  * (33) S. Khamseh, J.N. Roux, F. Chevoir, Phys. Rev. E 92(2), 022201 (2015)
  * (34) Y. Yamaguchi, S. Takada, T. Hatano, J. Phys. Soc. Jpn. 87(9), 094802 (2018)
  * (35) M. Badetti, A. Fall, F. Chevoir, J.N. Roux, Eur. Phys. J. E 41(5), 1 (2018)
  * (36) S. Mandal, M. Nicolas, O. Pouliquen, Phys. Rev. X 11(2), 021017 (2021)
  * (37) T.T. Vo, S. Nezamabadi, P. Mutabaruka, J.Y. Delenne, F. Radjai, Nat. Commun. 11(1), 1 (2020)
  * (38) T.T. Vo, P. Mutabaruka, S. Nezamabadi, J.Y. Delenne, F. Radjai, Phys. Rev. E 101(3), 032906 (2020)
  * (39) T.T. Vo, T. Nguyen-Thoi, Eur. Phys. J. E 43(10), 1 (2020)
  * (40) M. Macaulay, P. Rognon, Soft matter 17(1), 165 (2021)
  * (41) D.A. Head, Eur. Phys. J. E 22(2), 151 (2007)
  * (42) W. Zheng, H. Liu, N. Xu, Phys. Rev. E 94(6), 062608 (2016)
  * (43) A. Singh, V. Magnanimo, K. Saitoh, S. Luding, Phys. Rev. E 90(2), 022202 (2014)
  * (44) M. Macaulay, P. Rognon, J. Fluid Mech. 858 (2019)
  * (45) G. Lois, J. Blawzdziewicz, C.S. O’Hern, Phys. Rev. Lett. 100(2), 028001 (2008)
  * (46) D.J. Koeze, B.P. Tighe, Phys. Rev. Lett. 121(18), 188002 (2018)
  * (47) D.J. Koeze, L. Hong, A. Kumar, B.P. Tighe, Phys. Rev. Research 2, 032047 (2020)
  * (48) D. Evans, G. Morriss.  Non-equilibrium statistical mechanics of liquids (2008)
  * (49) S. Roy, A. Singh, S. Luding, T. Weinhart, Comput. Part. Mech. 3(4), 449 (2016)
  * (50) S. Roy, S. Luding, T. Weinhart, New J. Phys. 19(4), 043014 (2017)
  * (51) S. Roy, S. Luding, T. Weinhart, Phys. Rev. E 98(5), 052906 (2018)
  * (52) H. Shi, S. Roy, T. Weinhart, V. Magnanimo, S. Luding, Granul Matter. 22(1), 1 (2020)
  * (53) G. Lian, C. Thornton, M.J. Adams, J. Colloid Interface Sci. 161(1), 138 (1993)
  * (54) C.D. Willett, M.J. Adams, S.A. Johnson, J.P. Seville, Langmuir 16(24), 9396 (2000)
  * (55) M. Doi, S.F. Edwards, _The theory of polymer dynamics_ , vol. 73 (oxford university press, 1988)
  * (56) F. Da Cruz, S. Emam, M. Prochnow, J.N. Roux, F. Chevoir, Phys. Rev. E 72(2), 021309 (2005)
  * (57) J.D. Bernal, Proc. R. Soc. A 280(1382), 299 (1964)
  * (58) J.L. Finney, Proc. R. Soc. A 319(1539), 479 (1970)
  * (59) J.L. Finney, Proc. R. Soc. A 319(1539), 495 (1970)
  * (60) L. Oger, A. Gervois, J. Troadec, N. Rivier, Philos. mag., B 74(2), 177 (1996)
  * (61) R. Jullien, P. Jund, D. Caprion, D. Quitmann, Phys. Rev. E 54(6), 6035 (1996)
  * (62) R. Yang, R. Zou, A. Yu, Phys. Rev. E 65(4), 041302 (2002)
  * (63) T. Aste, T. Di Matteo, M. Saadatfar, T.J. Senden, M. Schröter, H.L. Swinney, EPL (Europhysics Letters) 79(2), 24003 (2007)
  * (64) J. Xu, R. Zou, A. Yu, Granular matter 9(6), 455 (2007)
  * (65) C. Rycroft, Voro++: A three-dimensional voronoi cell library in c++.  Tech. rep., Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States) (2009)