# MarkovGNN: Graph Neural Networks on Markov Diffusion

Md. Khaledur Rahman  Indiana University Bloomington<EMAIL_ADDRESS>,  Abhigya
Agrawal  Indiana University Bloomington<EMAIL_ADDRESS>and  Ariful Azad
Indiana University Bloomington<EMAIL_ADDRESS>
(2022)

###### Abstract.

Most real-world networks contain well-defined community structures where nodes
are densely connected internally within communities. To learn from these
networks, we develop MarkovGNN that captures the formation and evolution of
communities directly in different convolutional layers. Unlike most Graph
Neural Networks (GNNs) that consider a static graph at every layer, MarkovGNN
generates different stochastic matrices using a Markov process and then uses
these community-capturing matrices in different layers. MarkovGNN is a general
approach that could be used with most existing GNNs. We experimentally show
that MarkovGNN outperforms other GNNs for clustering, node classification, and
visualization tasks. The source code of MarkovGNN is publicly available at
https://github.com/HipGraph/MarkovGNN.

Graph Neural Network, Markov Clustering, Graph Clustering

††journalyear: 2022††copyright: acmcopyright††conference: Companion
Proceedings of the Web Conference 2022; April 25–29, 2022; Virtual Event,
Lyon, France††booktitle: Companion Proceedings of the Web Conference 2022 (WWW
’22 Companion), April 25–29, 2022, Virtual Event, Lyon, France††price:
15.00††doi: 10.1145/3487553.3524713††isbn: 978-1-4503-9130-6/22/04††ccs:
Computing methodologies Semi-supervised learning settings

##  1\. Introduction

Graph representation learning using graph neural networks (GNNs) is often
formulated by a message passing model (Gilmer et al., 2017). In this model, a
node $u$ receives messages from its direct neighbors and then updates $u$’s
embedding based on the received messages and its own features. While a GNN
with $L$ layers indirectly exchanges messages among $L$-hop neighbors,
messages in each layer are still restricted to direct neighbors. Message
passing via direct neighbors on the same static graph structure may limit the
expressive power of GNNs. Can we address this limitation by generating a
series of graphs (using a diffusion process) created from a given static
structure and then using the generated graphs at different layers of GNN? In
this paper, we answer this question affirmatively using a Markov diffusion
process (van Dongen, 2000) to generate a series of graphs. The goal is not to
develop a new GNN method, but to use Markov diffusion to improve the
performance of any existing GNN. When an existing GNN method uses different
Markov matrices at different layers, we call this augmented GNN a MarkovGNN.

Figure 1. Updating the embedding $x_{1}^{(2)}$ of $v_{1}$ at the second layer
of GNN. (a) For GCN, $x_{1}^{(2)}$ is computed directly from the embeddings of
$v_{2}$ and $v_{3}$ (shown by solid arrowheads), which indirectly depend on
the embeddings of $v_{1}$’s two-hop neighbors $v_{4}$ and $v_{5}$ (shown by
dashed arrowheads). (b) The second layer of MarkovGNN uses a different graph
structure with an edge $\\{v_{1},v_{4}\\}$ inserted and an edge
$\\{v_{3},v_{5}\\}$ deleted possibly because $v_{4}$ is in the same community
of $v_{1}$ whereas $v_{5}$ is in a different community. This graph structure
is updated based on a Markov process that captures the community pattern in
the graph. Thus, in MarkovGNN, $x_{1}^{(2)}$ is computed directly from the
embeddings of $v_{2}$, $v_{3}$, and $v_{3}$ (shown by solid arrowheads), which
are neighbors of $v_{1}$ in the current graph.

The central idea of MarkovGNN is to use a series of graphs at different layers
of a GNN so that different layers do not operate on the same static graph
structure. To understand the limitation of using the same static graph at
every layer, Fig. 1(a) considers the message passing at the second layer of a
graph convolutional network (GCN) (Kipf and Welling, 2017). At the second
layer of GCN, node $v_{1}$ receives information from 2-hop neighbors ($v_{4}$
and $v_{5}$) indirectly via 1-hop neighbors ($v_{2}$ and $v_{3}$). Such
indirect messages dilute information received from high-order neighbors and
may create the “over-smoothing” (Alon and Yahav, 2021; Oono and Suzuki, 2020)
problem resulting from hidden-layer representations being squeezed into fixed-
length vectors. What if we eliminate indirect messaging by using a new
structure created from the original graph? Fig. 1(b) shows a modified graph
where we added a new edge $\\{v_{1},v_{4}\\}$ and deleted an existing edge
$\\{v_{3},v_{5}\\}$ based on a diffusion model (here, to capture the community
structure in the original graph). If we use the modified graph shown in Fig.
1(b) at the second layer of a GNN, we promote direct messages (from $v_{4}$)
and suppress indirect messages (from $v_{5}$). We argue that using a series of
modified graphs in the hidden layers for direct messaging instead of a static
graph for indirect messaging benefits most existing GNN methods.

There are many ways to perturb a given graph to generate a series of graphs
such as the heat kernel and personalized PageRank. In this paper, we use a
modified Markov diffusion process (Van Dongen, 2000) because if its
flexibility in controlling diffusion while capturing the community structures
in the graph. The Markov process generates a sequence of graphs (called Markov
matrices) with an aim to find communities or clusters in the original graph (a
community represents a subset of nodes that are densely connected to each
other and loosely connected to other communities). Thus, the sequence of
Markov matrices promotes or demotes edges based on their relevance to the
underlying clustering pattern. Consequently, using Markov matrices at
different layers of a GNN helps learning from graphs that expose homophily
(Zhu et al., 2020) by organizing nodes into communities.

For graphs containing clustering patterns (e.g., graphs with high clustering
coefficients measured by counting the number of triangles in the graph),
MarkovGNN offers the following benefits. (1) By simulating the community
formation process, MarkovGNN directly captures community patterns into
different layers of a GNN. (2) Different Markov matrices increasingly put more
weight on intra-cluster edges and less weight on inter-cluster edges. Thus,
the Markov process (not the GNN) automatically controls information flows on
paths based on their relevance to form communities. This community-driven flow
eliminates the need to sample neighbors used in GraphSAGE (Hamilton et al.,
2017) or sample graphs used in GraphSAINT (Zeng et al., 2020). (3) As shown in
Fig. 1 (b), MarkovGNN only exchanges messages among 1-hop neighbors on
modified graphs. Hence, it can prevent the “over-smoothing” (Alon and Yahav,
2021; Oono and Suzuki, 2020) problem at high-degree nodes. (4) By pruning low
probable diffusion paths, the Markov process keeps the graph sparse at every
layer of the GNN and prevents the “neighbor-explosion” (Hamilton et al., 2017)
problem. (5) MarkovGNN is more expressive than Cluster-GCN (Chiang et al.,
2019) that precomputes clusters and uses the clustering at every layer of GCN.
We will show that a special variant of MarkovGNN is equivalent to Cluster-GCN.

This paper makes the following key contributions.

  1. (1)

We present MarkovGNN that uses a series of community-aware sparse adjacency
matrices at different layers of a GNN. MarkovGNN can augment most GNN methods,
and it is more general than GDC (Klicpera et al., 2019) and Cluster-GCN
(Chiang et al., 2019).

  2. (2)

By using residual connections, MarkovGNN supports deeper GNNs. MarkovGNN
reduces the neighbor-explosion problem in deeper layers by keeping Markov
matrices sparse.

  3. (3)

We conduct an extensive set of experiments to evaluate the performance of
MarkovGNN with respect to baseline GNN methods that use the same graph at
every layer. If vertices in the original graph are well-clustered, MarkovGNN
preserves the clustering pattern in the embedding space and improves the
performance of clustering, node classification, and visualization tasks.

##  2\. Background and Motivation

Let $G(V,E)$ denote a graph with $n=|V|$ vertices and $|E|$ edges. The
adjacency matrix is denoted by an $n\times n$ sparse matrix $\mathbf{A}$ where
$\mathbf{A}[i,j]=1$ if $\\{v_{i},v_{j}\\}{\in}E$, otherwise
$\mathbf{A}[i,j]=0$. Additionally, $\mathbf{X}\in\mathbb{R}^{n\times f}$
denotes the feature matrix where the $i$th row stores $f$-dimensional features
for the $i$th vertex. The $l$th convolution layer of a GNN computes embeddings
of every node based on the embeddings of neighbors as follows:

(1) $\mathbf{X}^{(l+1)}=\sigma(\mathbf{A}\mathbf{X}^{(l)}\mathbf{W}^{(l)}),$

where $\mathbf{X}^{(l)}{\in}\mathbb{R}^{n\times f_{l}}$ is $f_{l}$-dimensional
embeddings used as input to the $l$th layer,
$\mathbf{X}^{(l+1)}{\in}\mathbb{R}^{n\times f_{l+1}}$ is $f_{l+1}$-dimensional
embeddings generated as the output from the $l$th layer and
$\mathbf{W}^{(l)}{\in}\mathbb{R}^{f_{l}\times f_{l+1}}$ is the weight matrix.
In most instances of GNN, the adjacency matrix is normalized and ReLU is used
as the activation function.

The limitation of GNNs that use a fixed graph at every layer. If a static
graph is used at every layer of a GNN, it only passes messages between
adjacent nodes where the adjacency remains fixed throughout the execution of
the algorithm. While GCN and other message passing networks do leverage
higher-order neighborhoods in deeper layers, all messages must pass through
immediate neighbors which create bottlenecks when flowing messages from
higher-order neighbors (Alon and Yahav, 2021). To remove such bottlenecks,
high-pass and low-pass filters were proposed in a recent work with an aim to
control the information flow (Bo et al., 2021). However, in almost all GNNs,
the graph structure remains static at different layers. We argue that distinct
graph structures used in different layers of GNNs could benefit most existing
GNN methods. In particular, we aim to capture the formation of communities in
different layers of GNN to improve the expressivity of GNN models.

Previous work that uses diffusion in GNNs. The most notable approach is called
graph diffusion convolution (GDC) (Klicpera et al., 2019) where a diffusion
matrix $\tilde{\mathbf{A}}$ is constructed via a generalized graph diffusion
process:

(2)
$\tilde{\mathbf{A}}=\sum_{j=0}^{\infty}\beta_{j}\mathbf{T}^{j},\vspace{-0.17cm}$

where $\beta_{j}$ is the weighting coefficient and $\mathbf{T}$ is a column
stochastic adjacency matrix computed by $\mathbf{A}\mathbf{D}^{-1}$ and
$\mathbf{D}$ is the diagonal matrix of node degrees. The idea used in GDC is
interesting because of the well-known fact that graph diffusion can act as a
denoising filter similar to Gaussian filters on images (Egilmez et al., 2018).
However, GDC creates the diffusion matrix as a pre-processing step and uses
the same diffusion matrix $\tilde{\mathbf{A}}$ in every layer of a GNN. Hence,
GDC still uses the same connectivity matrix (albeit different from the
original graph) in every layer and thus misses the opportunity to
incrementally learn from the diffusion process. Furthermore, the diffusion
process typically results in a dense matrix, making the computation
prohibitively expensive. Even though GDC sparsifies $\tilde{\mathbf{A}}$
before using it in GNNs, the computation of $\mathbf{T}^{j}$ becomes the
bottleneck. Our goal is to capture the diffusion process (more specifically
via a Markov process) directly in different layers of GNN. We show that our
algorithm performs much better than GDC in identifying clustering patterns and
classifying nodes.

Related work. Graph representation learning is a well-studied research area
with numerous unsupervised and semi-supervised methods proposed in the
literature. Early unsupervised works focused on graph embedding based on
neighborhood exploration such as random-walk based methods (Perozzi et al.,
2014; Grover and Leskovec, 2016; Ribeiro et al., 2017; Rahman et al., 2020)
and $k$-hop neighborhood exploration-based methods (Tang et al., 2015). While
unsupervised embedding methods can capture the structure of a graph, they
often cannot utilize node features or domain-specific labels. Graph
Convolution Networks (GCN) is one of the first semi-supervised methods which
used graph convolutions for better representation learning (Kipf and Welling,
2017). Later, many other GNN methods have been proposed to tackle evolving
issues in GNNs and improve graph learning, e.g., GraphSAGE to incorporate
different types of aggregation strategies (Hamilton et al., 2017), Graph
Attention Networks (GAT) to consider the asymmetric influence of neighbors
(Veličković et al., 2018; Abu-El-Haija et al., 2018), FastGCN to improve the
training process (Chen et al., 2018), Cluster-GCN to resolve neighborhood
expansion problem and train deep networks (Chiang et al., 2019). Some recent
efforts pre-process the original graph and/or input features to create a
modified graph feeding it to all convolutional layers of the GNN (Klicpera et
al., 2019). We refer readers to a comprehensive survey for further details (Wu
et al., 2020). Most GNN methods use a static graph (often preprocessed) in all
GNN layers (Kipf and Welling, 2017; Hamilton et al., 2017; Klicpera et al.,
2019; Yang et al., 2020; Suresh et al., 2021) whereas some methods use a
mixing contribution of neighbors (Abu-El-Haija et al., 2019; Wu et al., 2019;
Liu et al., 2020). GraphSAINT uses different sampled subgraphs in different
minibatches but it still uses the same subgraph at every layer (Zeng et al.,
2020). Different from previous work, we put efforts to use different graph
structures in different convolutional layers.

##  3\. The Markov diffusion process

In this section, we review the Markov process (van Dongen, 2000) that
simulates flow within a graph such that the flow eventually gets trapped in
communities at the convergence of the process. The flow is simulated by random
walks with the assumption that higher-length paths between two nodes within a
community should stay within the community. The Markov process controls random
walks by computing transition probabilities between every pair of vertices.
The process incrementally promotes the probability of intra-cluster walks,
demotes the probability of inter-cluster walks, and removes paths with low
probabilities to discover natural clusters in the graph. Markov process
consists of four operations: (a) starting with a column (or row) stochastic
transition matrix and maintaining it until convergence, (b) expanding random
walks by the expansion operation, (c) strengthening intra-cluster walks and
weakening inter-cluster walks by the inflation operation, and (d) removes
paths with low probabilities by the prune operation. These operations are
repeated until converges.

The transition matrix. Let $\mathbf{D}$ be a diagonal matrix with
$\mathbf{D}[i,i]=\sum_{j=1}^{n}\mathbf{A}[i,j]$, where $\mathbf{A}$ is the
(weighted) adjacency matrix of the graph. The Markov process starts with a
column-stochastic matrix $\mathbf{M}=\mathbf{A}\mathbf{D}^{-1}$. The $i$th
column of $\mathbf{M}$ sums to one, and $\mathbf{M}[i,j]$ is interpreted as
the probability of following the edge from vertex $j$ to vertex $i$ in the
random walk.

The expansion operation. The Markov process expands random walks by squaring
the Markov matrix $\mathbf{M}$. Thus, the expansion operation in the $i$th
iteration computes $\mathbf{M_{i-1}}{\times}\mathbf{M_{i-1}}$ where
$\mathbf{M_{i-1}}$ is the output of the previous iteration. The expansion step
is similar to the generalized graph diffusion (Klicpera et al., 2019) and
typically converges to an idempotent matrix according to the Perron–Frobenius
theorem.

The inflation operation. The parameterized inflation operation with parameter
$r$ performs Hadamard power of the matrix denoted by $\mathbf{M}^{\circ r}$
and makes the matrix column stochastic again as follows:

(3) $\mathbf{M}^{\circ
r}[i,j]=\mathbf{M}[i,j]^{r}/\sum_{i=1}^{n}\mathbf{M}[i,j]^{r}$

The goal of the non-linear inflation step is to strengthen higher probability
paths (intra-cluster connections) and weaken lower probability paths (inter-
cluster connections).

The prune operation. To keep the Markov matrices sparse, we prune small
entries that are below a threshold $\theta$. This step ensures that random-
paths are contained within communities and that the matrix does not become
dense.

Algorithm 1 Generating Markov matrices

Input: $\mathbf{A}$: adjacency matrix of the graph, $r:$ value of element-wise
inflation, $\theta:$ pruning threshold

Output: $\mathbf{M}_{1},\mathbf{M}_{2},...,\mathbf{M}_{k}$: column stochastic
matrices from all steps of the Markov process.

1:procedure Markov($\mathbf{A},r,\theta$)

2: $\mathbf{M_{1}}\leftarrow\text{ColStochastic}(\mathbf{A})$

3: for $k\leftarrow 2..$ do

4: $\mathbf{M}_{k}\leftarrow\mathbf{M}_{k-1}\times\mathbf{M}_{k-1}$
$\triangleright$ the expansion operation

5: $\mathbf{M}_{k}\leftarrow\text{ColStochastic}(\mathbf{M}_{k}^{\circ
r})$$\triangleright$ inflation

6:
$\mathbf{M}_{k}\leftarrow\text{ColStochastic}(\text{Prune}(\mathbf{M}_{k},\theta))$
$\triangleright$ Sparsify the matrix by pruning small entries less than
$\theta$

7: if $\mathbf{M}_{k}\approx\mathbf{M}_{k-1}$ then

8: return $\mathbf{M}_{1},\mathbf{M}_{2},...,\mathbf{M}_{k}$

Generating Markov matrices for GNNs. Algorithm 1 shows the steps used to
generate Markov diffusion matrices. The algorithm starts with a column
stochastic matrix $\mathbf{M}_{1}$ where each column sums to 1 and maintains
column stochastic matrices in each iteration. The $k$th iteration performs the
expansion, inflation, and pruning operations. The algorithm converges when the
stochastic matrix does not change in successive iterations. We store all
stochastic matrices $\mathbf{M}_{1},\mathbf{M}_{2},...,\mathbf{M}_{k}$ from
different iterations so that they can be used in different layers of a GNN.
While there is a memory overhead to store $k$ Markov matrices, it makes the
GNN training much faster as we will show in the next section.

Figure 2. The bottom panel shows column stochastic matrices generated by the
Markov process for hypothetical graphs with bi-directional edges shown by
lines. Here, $\mathbf{M}_{1}$ is the input stochastic matrix representing the
given graph, $\mathbf{M}_{5}$ is the converged matrix with four discovered
communities, and $\mathbf{M}_{2}$, $\mathbf{M}_{3}$, and $\mathbf{M}_{4}$ are
intermediate matrices capturing the formation of communities. As communities
are discovered, edges are added and deleted in different iterations. The top
panel shows a GNN that uses a subset of Markov matrices in three convolutional
layers where $\alpha$ is the contribution factor of residual connections.

Convergence. It has already been shown that the Markov process converges
quadratically (Van Dongen, 2000). In most practical graphs, the process
converges in less than 30 iterations. Upon convergence, connected components
are treated as clusters in the original graph.

Hyperparameters in the Markov process. The exponent $r$ in inflation and the
pruning threshold $\theta$ are hyperparameters that could be tuned to obtained
different matrices. A higher value of the parameter $r$ demotes flow along
long paths, which creates fine-grained communities. Thus, the number and
granularity of final communities can be controlled by $r$. When $r$ is set to
1 (i.e., no inflation), the Markov process is equivalent to a general
diffusion process (Klicpera et al., 2019). $\theta$ controls the sparsity of
Markov matrices. We empirically identify suitable values of these parameters.

##  4\. The MarkovGNN Algorithm

The MarkovGNN algorithm exchanges the static adjacency matrix with a series of
Markov matrices at different layers of a GNN. This approach could be used with
most popular GNN methods such as GCN, GAT, and GraphSAINT. Figure 2
illustrates the forward propagation of a MarkovGNN. In this example, the GNN
has three layers (shown in the top panel), but the Markov process converges in
four iterations (in the bottom panel). Hence, layer 1, 2, and 3 of GNN use
$\mathbf{M}_{1}$, $\mathbf{M}_{3}$, and $\mathbf{M}_{5}$, respectively.

Algorithm 2 MarkovGCN forward propagation

Input: $\mathbf{A}$: adjacency matrix of an input graph, $\mathbf{X}$: input
feature, $L:$ number of layers, $r$: inflation parameter, $\theta$: pruning
threshold

Output: $\mathbf{X}^{(L+1)}$: node representation from the last layer.

1:procedure MarkovGCN($\mathbf{A},L,r,\theta$)

2:
$\\{\mathbf{M}_{1},\mathbf{M}_{2},...,\mathbf{M}_{k}\\}\leftarrow\text{MARKOV}(\mathbf{A},r,\theta)$
$\triangleright$ Alg. 1

3: $\mathbf{X}^{(1)}\leftarrow\mathbf{X}$

4: $\mathbf{W}^{(1)},...\mathbf{W}^{(L)}$ are weight matrices

5: $\\{\mathbf{A}^{(1)},\mathbf{A}^{(2)},...,\mathbf{A}^{(L)}\\}\leftarrow$
select $L$ matrices from
$\\{\mathbf{M}_{1},\mathbf{M}_{2},...,\mathbf{M}_{k}\\}$ one for each layer
$\triangleright$ assuming $L{<=}k$

6: for $l=1$ to $L$ do

7:
$\mathbf{X}^{(l+1)}=\alpha\sigma(\mathbf{A}^{(l)}\mathbf{X}^{(l)}\mathbf{W}^{(l)})+(1-\alpha)\phi(\mathbf{X}_{1}|\ldots|\mathbf{X}_{l-1})$

8: return $\mathbf{X}^{(L+1)}$

Residual connections. A MarkovGNN needs more layers than a typical GNN to
utilize various Markov matrices at different layers. To train such deep GNN
models, we use residual connection (He et al., 2016) which has also been used
in other deep networks (Li et al., 2021; Chen et al., 2020). With residual
connections, the embedding in the $l$-th layer of MarkovGNN is computed as
follows:

(4)
$\mathbf{X}^{(l+1)}=\alpha\sigma(\mathbf{A}^{(l)}\mathbf{X}^{(l)}\mathbf{W}^{(l)})+(1-\alpha)\phi(\mathbf{X}_{1}|\ldots|\mathbf{X}_{l-1})$

where, $\sigma$ is the activation function, $\phi$ is a selector and $\alpha$
is a hyper-parameter that controls the contribution of the residual
connections. Note that the operator $\phi$ in $l$-th layer selects only one
representation matrix from the previous $l-1$ layers.

The MarkovGNN algorithm. Algorithm 2 describes the forward propagation of the
MarkovGCN algorithm. If a GNN has $L$ layers, we select $L$ matrices
$\\{\mathbf{A}^{(1)},\mathbf{A}^{(2)},...,\mathbf{A}^{(L)}\\}$ from the
collection of stochastic matrices returned by Algorithm 1. Then, in the $l$th
layer of the GNN, convolution is performed on $\mathbf{A}^{(l)}$ instead of
the original adjacency matrix (line 7 of Algorithm 2). Here, we assume that
the number of iterations $k$ needed for the convergence of the Markov process
is greater than $L$. This is a reasonable assumption since $k$ is typically
more than 20 for most graphs we considered, and most GNNs have less than 10
layers in practice. The selection of $L$ matrices from $k$ available matrices
is an algorithmic choice for which we show various experimental results.
However, we always present Markov matrices to GNN layers in the same order
they are generated in the Markov process.

Training MarkovGNN. In Algorithm 2, $\sigma$ represents an activation function
in line 7 for which either ReLU or LeakyReLU is a common choice. We add a
tunable dropout parameter to each convolutional layer. In addition, a softmax
layer is added to the last layer of MarkovGNN. Similar to other GNNs, our goal
is to learn weight matrices by optimizing a loss function as follows:

(5) $\mathcal{L}=\sum_{i\in\mathrm{Y}_{L}}loss(y_{i},\mathbf{X}^{(L)}_{i}),$

where, $\mathrm{Y}_{L}$ represents all vertices that have ground truth labels,
and $\mathbf{X}^{(L)}_{i}$ represents $i$th row of $\mathbf{X}^{L}$ having
ground truth label $y_{i}$. In MarkovGNN, we use the standard negative log-
likelihood loss function and the Adam optimizer (Kingma and Ba, 2015) to train
models.

###  4.1. Variants of MarkovGNN

MarkovGNN is a general technique that could be used with any GNN. For example,
when Markov matrices are used with GCN, we call this combination MarkovGCN
(similarly, MarkovGAT, MarkovGraphSAINT, etc.). Furthermore, the Markov
matrices can be combined differently to obtain several interesting variants of
GNNs as shown in Table 1 and described below.

Table 1. GNN variants by combining Markov matrices. $\mathbf{A}^{(1)}$ |  $\mathbf{A}^{(2)}$… | $\mathbf{A}^{(L)}$ | Equivalent GNN models
---|---|---|---
$\mathbf{M}_{1}$ | $\mathbf{M}_{1}$ | $\mathbf{M}_{1}$ | GCN (Kipf and Welling, 2017)
$\sum_{i=1}^{k}\mathbf{M}_{i}$ | $\sum_{i=1}^{k}\mathbf{M}_{i}$ | $\sum_{i=1}^{k}\mathbf{M}_{i}$ | GDC (Klicpera et al., 2019)
$\mathbf{M}_{k}$ | $\mathbf{M}_{k}$ | $\mathbf{M}_{k}$ | Cluster-GCN (Chiang et al., 2019)
$\mathbf{M}_{1}$ | $\mathbf{M}_{2}$ | $\mathbf{M}_{k}$ | General Markov-GNN

uMarkovGNN: union of Markov matrices. In row 2 of Table 1 we use
$\mathbf{\tilde{M}}{=}\sum_{i=1}^{k}\mathbf{M}_{i}$ as a unified Markov matrix
and use it in all layers of a GNN. For example, the forward propagation at
layer $l$ of GCN could be computed by
$\mathbf{\tilde{M}}\mathbf{X}^{(l)}\mathbf{W}^{(l)}$ where the graph structure
$\mathbf{\tilde{M}}$ does not change in different layers. We call this
approach uMarkovGNN which is equivalent to the diffusion matrix considered in
Eq. 2. However, uMarkovGCN may perform better than GDC because the former
preserves a unified view of communities.

MarkovClusterGNN. In row 3 of Table 1 we use the converged Markov matrix
$\mathbf{M}_{k}$ at every layer of GNN (that is, the computation at layer $l$
is $\mathbf{M}_{k}\mathbf{X}^{(l)}\mathbf{W}^{(l)}$). We call this approach
MarkovClusterGNN, which is equivalent to Cluster-GCN (Chiang et al., 2019).
Cluster-GCN also identifies clusters using a graph partitioner and then uses
the clustered graph at every layer of GCN.

###  4.2. Properties of MarkovGNN

Computational complexity. Let $\mathbf{M}_{k}$ have $d_{k}$ nonzeros per
column (that is, the average degree of the corresponding graph is $d_{k}$).
Then, the expected cost of computing $\mathbf{M}_{k}$ in line 4 of Algorithm 1
is $nd_{k}^{2}$ (Ballard et al., 2013). Similarly, let $\mathbf{A}^{(l)}$ have
$d_{l}$ nonzeros per column on average. Then the computational cost of the
forward propagation is

(6) $\sum_{l=1}^{L}(nd_{l}f_{l}+nf_{l}f_{l+1}).\vspace{-6pt}$

The role of sparsity. As communities are formed in the Markov process, the
corresponding matrices become increasingly sparse. As a result, $d_{l+1}$ is
typically smaller than $d_{l}$ in deeper layers of MarkovGNN. Thus, MarkovGNN
is expected to run faster than other other diffusion based methods such as
GDC.

An upper bound on the number of layers $L$. In most GNN methods, the number of
layers is selected in an ad hoc fashion. By contrast, MarkovGNN sets an upper
bound on the number of layers to the number of Markov matrices ($k$) returned
by Algorithm 1.

A potential solution to the neighbor-explosion problem. As communities are
formed, the Markov process removes intra-cluster edges via the inflation and
pruning operations. These steps keep the graph sparse at every layer and could
prevent the neighbor-explosion problem observed in GNNs.

The role of sampling from
$\\{\mathbf{M}_{1},\mathbf{M}_{2},...,\mathbf{M}_{k}\\}$ in GNN layers. The
selection of $L$ matrices from $k$ available Markov matrices is an algorithmic
choice for which we can consider various selection strategies. This selection
of matrices from a pool of available matrices makes it unnecessary to sample
neighbors used in GraphSAGE (Hamilton et al., 2017) or to sample minibatches
used in GraphSAINT (Zeng et al., 2020).

Limitations of MarkovGNN. Since the Markov process iteratively identifies the
community structure in the graph, MarkovGNN is expected to excel when the
graph has well-defined communities. Our experimental results show clear
evidence in favor of this hypothesis. However, when a graph has no clear
clustering patterns MarkovGNN may not perform better than other GNN methods.

##  5\. Experiments

###  5.1. Experimental Setup

Overview of experiments. We conduct an extensive set of experiments by
coupling Markov matrices with GCN, GAT, and GraphSAINT. As MarkovGNN is
expected to capture the clustering patterns in a graph, we compare the
clustering effectiveness of MarkovGNN with other GNN and unsupervised
algorithms. We also examine the effectiveness of MarkovGNN in node
classification and graph visualization tasks.

Experimental platforms. We conduct all experiments on an Intel Skylake
processor with 48 cores and 128GB memory. We have implemented MarkovGNN using
PyTorch Geometric (PyG) (Fey and Lenssen, 2019) version 1.6.1 built on top of
PyTorch 1.7.0. We used GCN, GAT, Cluster-GCN, GDC, and GraphSAINT
implementations from PyG. We run all experiments 10 times and report the
average results.

Table 2. Summary of Datasets. ‘-’ indicates the absence of input features. Graphs are sorted based on the descending order of the Average Clustering Coefficient (ACC). Graphs | Nodes | Edges | Classes | Features | ACC
---|---|---|---|---|---
USAir | 1,190 | 13,599 | 4 | - | 0.50
Chameleon | 2,277 | 36,101 | 3 | 3,132 | 0.48
Squirrel | 5,201 | 217,073 | 3 | 31,48 | 0.42
Email | 1,005 | 25,571 | 42 | - | 0.40
Cora | 2,708 | 10,556 | 7 | 1433 | 0.24
Citeseer | 3,327 | 9,104 | 6 | 3703 | 0.14
PPI | 2,361 | 7,182 | 13 | - | 0.13
Actor | 7,600 | 33,544 | 5 | 931 | 0.08

Datasets. We experimented with eight graphs from different domains such as
actor co-occurrence network in films (Actor), airport networks (USAir),
community network (Email), biological network (PPI), wikipedia networks
(Chameleon and Squirrel), and citation networks (Cora and Citeseer). All these
graphs have been used previously in different research papers (Kipf and
Welling, 2017; Rozemberczki et al., 2021; Ribeiro et al., 2017; Pei et al.,
2020; Bo et al., 2021). Table 2 provides a summary of the graphs. The
clustering coefficient of a node in a network measures the propensity to form
a cluster and the Average Clustering Coefficient (ACC) is a generalization
that represents this score for the whole network (Saramäki et al., 2007). For
networks without input features, we use a one-hot encoding representation of
vertices and feed it to the GNNs as input features.

Experiment settings. Among the graphs having input features, Cora and Citeseer
have high assortativity whereas Chameleon and Squirrel have low assortativity
(Bo et al., 2021). Thus, for the classification task using Cora and Citeseer
networks, we use 30 vertices per class for training, 500 vertices for
validation, and 1000 vertices for testing. For all other networks, we use 70%
vertices for training, 10% vertices for validation, and the remaining 20%
vertices for testing.

In clustering experiments, we infer the predicted labels of all nodes using
the trained models and compute clustering measures such as the Adjusted Rand
Index. Similarly, for visualization, we infer 64-dimensional embeddings from
the last convolutional layer of GNN models and generate 2D layouts using t-SNE
(Van der Maaten and Hinton, 2008). Table 5 reports a list of parameters used
to run GNN methods.

Figure 3. Comparing the quality of predicted clustering using the ARI index
(higher is better). In each plot, bars are organized in three groups based on
their base GNN models: GCN (first four bars), GAT (5th and 6th bars), and
GraphSAINT (the last two bars). MarkovGNN represents the last bar in each
group, and it generally improves the performance of the base GNN model.

Baselines. We choose five semi-supervised methods (GCN, GDC, GAT, Cluster-GCN,
and GraphSAINT) to compare the effectiveness of proposed method. GCN is the
pioneering work using graph convolutions in GNN models (Kipf and Welling,
2017). GAT uses an attention mechanism by considering the influence of
neighbors of a vertex in a graph (Veličković et al., 2018). GDC uses diffusion
on a graph (as shown in Eq. 2) and then feeds the modified graph to GNNs
(Klicpera et al., 2019). In our analysis, we consider GCN as an underlying
method for GDC. Cluster-GCN partitions the graph using a clustering technique
before feeding it to GNN (Chiang et al., 2019). GraphSAINT is a recent method
based on different types of sampling approaches such as random-walk sampling
and edge sampling (Zeng et al., 2020). Besides these semi-supervised methods,
we also experimented with two unsupervised graph embedding methods such as
DeepWalk (Perozzi et al., 2014) and struc2vec (Ribeiro et al., 2017). We set
different hyper-parameters of other GNN models in PyG framework. Tables 5 and
6 provide a summary of the parameters used in our experiments.

###  5.2. Clustering Performance of MarkovGNN

Theoretically, we expect that MarkovGNN would detect the clustering pattern
present in the original graphs. To validate this expectation, we use predicted
vertex labels from MarkovGNN and other GNN methods to partition the graph into
predicted clusters. We then use true class labels as ground truth clusters and
measure the similarity between predicted and ground-truth clustering patterns.

Cluster similarity metrics. To measure the similarity between the predicted
and given clusterings, we use the adjusted Rand index (ARI) (Hubert and
Arabie, 1985) and V-Measure (Rosenberg and Hirschberg, 2007). The original
Rand index identifies (1) the number of vertex pairs that are in the same
cluster in both given and predicted clusterings (true positives) and (2) the
number of vertex pairs that are in different clusters in given and predicted
clusterings (true negatives). The true positives and true negatives are added
and normalized by the total number of vertex pairs to get the Rand index. ARI
simply adjusts the Rand index for the chance grouping of elements. The higher
the ARI the more similar two clusterings are. The V-Measure (Rosenberg and
Hirschberg, 2007) is computed from the homogeneity ($h$) and completeness
($c$) scores. Homogeneity computes how different the predicted labels are for
each true cluster and completeness computes how different the clusters are for
each type of predicted label. The V-Measure computes the weighted harmonic
mean of $h$ and $c$ i.e., $\frac{(1+\beta)hc}{\beta h+c}$, where $\beta{>}1$
puts more weight on homogeneity and $\beta{<}1$ puts more weight on
completeness (in our case, $\beta{=}1$). A higher value of V-Measure
represents more similar clusterings.

Figure 4. Comparing the performance in classifying nodes using different
variants of MarkovGNN and other GNN methods. In each plot, bars are organized
in three groups based on their base GNN models: GCN (first four bars), GAT
(5th and 6th bars), and GraphSAINT (the last two bars). MarkovGNN represents
the last bar in each group, and it often improves the performance of the base
GNN model.

MarkovGNN improves the prediction of clusters. Fig. 3 compares the quality of
clustering predictions (using ARI) from various GNN methods. We observe that
the Markov process improves the clustering performance of GCN, GAT, and
GraphSAINT for almost all graphs. For example, MarkovGCN performs better than
GCN for all graphs, MarkovGAT performs better than GAT for all graphs except
Cora and Citeseer, and MarkovGraphSAINT performs better than GraphSAINT for
all graphs except Cora and Citeseer. Furthermore, MarkovGCN consistently
outperforms GDC and Cluster-GCN even though the latter two methods use
diffusion and clustering to preprocess the graph. Overall, for any graph,
there is at least one variant of MarkovGNN that performs better than any other
approach according to the ARI measure. For example, MarkovGCN performs the
best for Cora even though Markov matrices do not improve the performance of
GraphSAINT and GAT on this graph. MarkovGNN also exhibits better clustering
performance according to V-Measure (expect for the PPI network where
GraphSAINT performs the best) as shown in Table 3. We also experimented with
unsupervised embedding methods DeepWalk, and struc2vec for community detection
(see Fig. 8). However, the unsupervised methods perform much worse than GNNs,
which is not surprising because unsupervised methods do not use class labels
when finding the embedding or communities.

Table 3. V-Measure (avg. of 10 runs) for different methods. A higher value of V-Measure means a better prediction of clusters. The best value for each dataset is highlighted in green. Clust.GCN - Cluster-GCN.  Methods |  USAir |  Cham. |  Squi. |  Email |  Cora |  Cite. |  PPI |  Actor
---|---|---|---|---|---|---|---|---
GCN |  0.564 |  0.436 |  0.297 |  0.778 |  0.616 |  0.418 |  0.534 |  0.211
GAT |  0.443 |  0.479 |  0.152 |  0.653 |  0.636 |  0.425 |  0.318 |  0.014
Clust.GCN |  0.591 |  0.532 |  0.38 |  0.86 |  0.588 |  0.383 |  0.609 |  0.058
GDC |  0.263 |  0.414 |  0.26 |  0.634 |  0.629 |  0.418 |  0.359 |  0.091
GraphSAINT |  0.448 |  0.566 |  0.217 |  0.817 |  0.517 |  0.315 |  0.673 |  0.221
MarkovGCN |  0.607 |  0.576 |  0.493 |  0.86 |  0.638 |  0.426 |  0.601 |  0.228

When does MarkovGNN excel in predicting clusters? When a graph has a high
clustering coefficient, MarkovGNN predicts the original clustering much better
than its peers. This behavior can be observed in Fig. 3 and Table 3. For
example, in Fig. 3, MarkovGCN performs much better than GCN for the USAir,
Chameleon, Squirrel, and Email graphs that have clustering coefficients
greater than or equal to $0.4$ (see Table 2). Cluster-GCN also performs well
(slightly worse than MarkovGCN) for these graphs, indicating the usefulness of
cluster-based GNN methods. However, for Cora, Citeseer, and PPI that have a
clustering coefficient less than 0.25, MarkovGCN performs slightly better than
GCN. Hence, the utility of MarkovGNN may diminish when the input graph lacks
clustering patterns. MarkovGNN not only attains higher values for ARI and
V-measure, it also achieves higher true positive rates for individual classes.
See Table 7 in the appendix for an example.

Figure 5. 2D layouts of the Chameleon graph from the embedding obtained from
GCN in the left and MarkovGCN in the right. 2D layouts are obtained using
t-SNE.

Figure 6. Parameter sensitivity of MarkovGCN showing (a) influence of pruning
threshold on Citeseer network with corresponding ARI and the remaining edges
for the 2nd convolutional layer, (b) influence of inflation parameter for
Chameleon and PPI networks, and (c) effect of using row vs. column stochastic
matrices in Algorithm 1.

###  5.3. Node Classification

Fig. 4 shows the test accuracy in predicting node labels from different GNN
methods. We observe that at least one variant of MarkovGNN shows the best test
accuracy for all graphs except for Chameleon and Squirrel where GraphSAINT and
Cluster-GCN, respectively, perform better than any other methods. Generally,
MarkovGNN improves the performance of its base GNN model in most graphs. In
particular, MarkovGCN tends to classify nodes well for graphs where graph-
structural properties significantly impact node classification such as USAir.
Similarly, the Email network has known ground truth communities that MarkovGCN
can exploit to outperform other methods significantly. For widely studied
graphs such as Cora and Citeseer, MarkovGCN marginally outperforms GCN.
However, we do not expect significant improvement for Cora and Citeseer
because of their lower average clustering coefficient values as shown in Table
2. Note that GraphSAINT is a recent method that performs exceptionally well
for Chameleon and Squirrel networks. However, GraphSAINT does not perform well
(often worse than vanilla GCN) for graphs with known communities (e.g., the
Email network) and structures (e.g., USAir). By contrast, these are the
networks where MarkovGCN excels.

###  5.4. Visualization of clusters

If an embedding of a graph captures the community structure, it should provide
an aesthetically pleasing visualization by keeping nodes well clustered in the
embedding space. To show the impact of MarkovGNN, we take the 64-dimensional
embedding from the last layer of a GNN and use t-SNE (Van der Maaten and
Hinton, 2008) to project 64-dimensional embeddings to 2D space. We show the
visualizations for the Chameleon graph from GCN and MarkovGCN in Fig. 5, where
nodes from the same ground-truth class are shown with the same color. Fig. 5
shows that the visualization obtained from MarkovGCN tends to place vertices
from the same class together while keeping different classes well separated.
For example, we can easily identify three clusters from the MarkovGCN-based
visualization in the right subfigure, but the visualization from GCN failed to
separate clusters from one another.

Table 4. Average runtimes of 10 runs in seconds are reported. For all GNNs, we only measure the training time of 200 epochs in the PyG framework.  Graphs |  GCN |  GAT |  GDC |  Gr.SAINT |  MarkovGCN
---|---|---|---|---|---
USAir |  7.4 |  21.5 |  95.3 |  212.4 |  7.7
Email |  14.3 |  31.8 |  98.2 |  251.9 |  11.3
PPI |  13.7 |  32.6 |  228.6 |  167.2 |  17.3
Chamel. |  26.8 |  60.3 |  222.3 |  441.1 |  30.6
Squirrel |  229.9 |  327.3 |  678.9 |  2,110.9 |  188.8
Cora |  10.9 |  31.2 |  245.8 |  160.8 |  16.4
Citeseer |  24.1 |  67.9 |  317.8 |  220.3 |  41.7

###  5.5. Running Time

We measure the average training time of all methods for 10 runs (available in
Table 4). Despite needing additional time in generating Markov matrices
(Algorithm 1), MarkovGCN and GCN have comparable running time. On the other
hand, GDC runs a diffusion step that is computationally similar to Markov
clustering. However, effective sparsification by inflation and pruning makes
MarkovGCN an order of magnitude faster than GDC. Similarly, MarkovGCN runs
$10\times$ to $100\times$ faster than GraphSAINT. The advantage of MarkovGCN
is originated from the decreasing number of edges in deeper layers. By
contrast, other GNN methods use the same graph in all layers.

###  5.6. Parameter Sensitivity

The impact of the pruning threshold ($\theta$). The pruning threshold plays an
important role in the convergence of the Markov process. Fig. 6 (a) reports
the ARI of MarkovGCN using Citeseer for different values of $\theta$. For the
Citeseer network, MarkovGCN achieves the best ARI score when $\theta\approx
0.15$. We also observe in Fig. 6 (a) that an increasing threshold decreases
the number of edges in the second layer of MarkovGCN. Thus, the computation
can be made much faster by sacrificing the accuracy slightly. If ARI or other
measures do not change for different values of $\theta$, we prefer to use
higher thresholds which typically make training and inference faster.

Impact of the inflation parameter ($r$). To find an optimal value of inflation
in Algorithm 1, we vary $r$ from 1.1 to 3.1 and run the experiments using
Chameleon and PPI networks. Fig. 6(b) shows that MarkovGCN is less sensitive
to the inflation parameter. We observe better clustering results when $r$ is
close to 1.5 and keep this value as a default value of $r$.

Column vs row stochastic matrices. Algorithm 1 maintains column stochastic
matrices in every iteration. Fig. 6(c) shows that row stochastic matrices also
perform similarly for different graphs.

Figure 7. Layer-wise validation accuracy in GCN model for the Cora network.

Residual connections with a hyper-parameter ($\alpha$). In Eqn. 4, we aim to
control the residual connections with a hyper-parameter $\alpha$. To select a
suitable residual connection for $\phi$ in Eqn. 4, we conduct an experiment to
predict the impact of different layers on the predictive performance of a
four-layer GNN. Following the protocol used by Raghu et al. (Raghu et al.,
2017), we disabled gradient updates for all layers except one and plot the
observed validation accuracy for different layers in Fig. 7. For the Cora
network, we observe that the input layer (layer 1) has the highest accuracy
followed by layer 4. The contribution of layer 1 is expected to be significant
for graphs such as Cora that include node features. Thus, when adding residual
connections, we select layer 1 in Eqn. 4.

##  6\. Conclusions

Communities in a network can be discovered by probabilistic random walks with
pruning that iteratively promote intra-cluster edges and demote inter-cluster
edges. The proposed algorithm MarkovGNN captures snapshots of this community
discovery process in different layers of graph neural networks. We conclude
with experimental evidence that using Markov matrices in different layers
improve the performance of clustering, node classification and visualization
tasks. MarkovGNN excels on graphs having well-defined community structures
(e.g., a high value of the clustering coefficient) because it can preserve the
communities in the embedding space. However, if a graph lacks a community
structure, MarkovGNN performs no worse than GCN. Overall, MarkovGNN brings in
the philosophy of using different community-capturing graphs at different
layers of a GNN and provides a flexible design space for graph representation
learning.

###### Acknowledgements.

This research was supported by the Applied Mathematics Program of the DOE
Office of Advanced Scientific Computing Research under contract number DE-
SC0022098.

## References

  * (1)
  * Abu-El-Haija et al. (2018) Sami Abu-El-Haija, Bryan Perozzi, Rami Al-Rfou, and Alex Alemi. 2018.  Watch your step: Learning node embeddings via graph attention.  _Proceedings of NeurIPS_ (2018).
  * Abu-El-Haija et al. (2019) Sami Abu-El-Haija, Bryan Perozzi, Amol Kapoor, Nazanin Alipourfard, Kristina Lerman, Hrayr Harutyunyan, Greg Ver Steeg, and Aram Galstyan. 2019.  Mixhop: Higher-order graph convolutional architectures via sparsified neighborhood mixing. In _international conference on machine learning_. PMLR, 21–29.
  * Alon and Yahav (2021) Uri Alon and Eran Yahav. 2021.  On the Bottleneck of Graph Neural Networks and its Practical Implications. In _Proceedings of ICLR_.
  * Ballard et al. (2013) Grey Ballard, Aydin Buluc, James Demmel, Laura Grigori, Benjamin Lipshitz, Oded Schwartz, and Sivan Toledo. 2013.  Communication optimal parallel multiplication of sparse random matrices. In _Proceedings of the Twenty-fifth annual ACM Symposium on Parallelism in Algorithms and Architectures_. ACM, 222–231.
  * Bo et al. (2021) Deyu Bo, Xiao Wang, Chuan Shi, and Huawei Shen. 2021\.  Beyond Low-frequency Information in Graph Convolutional Networks. In _Proceedings of AAAI_.
  * Chen et al. (2018) Jie Chen, Tengfei Ma, and Cao Xiao. 2018.  FastGCN: Fast learning with graph convolutional networks via importance sampling.  _Proceedings of ICLR_ (2018).
  * Chen et al. (2020) Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. 2020.  Simple and deep graph convolutional networks. In _International Conference on Machine Learning_. PMLR, 1725–1735.
  * Chiang et al. (2019) Wei-Lin Chiang, Xuanqing Liu, Si Si, Yang Li, Samy Bengio, and Cho-Jui Hsieh. 2019\.  Cluster-gcn: An efficient algorithm for training deep and large graph convolutional networks. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 257–266.
  * Egilmez et al. (2018) Hilmi E Egilmez, Eduardo Pavez, and Antonio Ortega. 2018\.  Graph learning from filtered signals: Graph system and diffusion kernel identification.  _IEEE Trans. on Signal and Information Processing over Networks_ 5, 2 (2018), 360–374.
  * Fey and Lenssen (2019) Matthias Fey and Jan Eric Lenssen. 2019.  Fast graph representation learning with PyTorch Geometric.  _arXiv preprint arXiv:1903.02428_ (2019).
  * Gilmer et al. (2017) Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. 2017\.  Neural message passing for quantum chemistry. In _International conference on machine learning_. PMLR, 1263–1272.
  * Grover and Leskovec (2016) Aditya Grover and Jure Leskovec. 2016.  node2vec: Scalable feature learning for networks. In _Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining_. 855–864.
  * Hamilton et al. (2017) William L Hamilton, Rex Ying, and Jure Leskovec. 2017\.  Inductive representation learning on large graphs. In _Proceedings of the 31st International Conference on Neural Information Processing Systems_. 1025–1035.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2016\.  Deep residual learning for image recognition. In _Proceedings of the IEEE conference on computer vision and pattern recognition_. 770–778.
  * Hubert and Arabie (1985) Lawrence Hubert and Phipps Arabie. 1985.  Comparing partitions.  _Journal of classification_ 2, 1 (1985), 193–218.
  * Kingma and Ba (2015) Diederik P Kingma and Jimmy Ba. 2015.  Adam: A method for stochastic optimization. In _Proceedings of ICLR_.
  * Kipf and Welling (2017) Thomas N Kipf and Max Welling. 2017.  Semi-supervised classification with graph convolutional networks. In _Proceedings of ICLR_.
  * Klicpera et al. (2019) Johannes Klicpera, Stefan Weißenberger, and Stephan Günnemann. 2019\.  Diffusion improves graph learning.  _Advances in Neural Information Processing Systems_ 32 (2019), 13354–13366.
  * Li et al. (2021) Bing Li, Zhi Li, and Yilong Yang. 2021.  Residual attention graph convolutional network for web services classification.  _Neurocomputing_ 440 (2021), 45–57.
  * Liu et al. (2020) Meng Liu, Hongyang Gao, and Shuiwang Ji. 2020.  Towards deeper graph neural networks. In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 338–348.
  * Oono and Suzuki (2020) Kenta Oono and Taiji Suzuki. 2020.  Graph neural networks exponentially lose expressive power for node classification. In _Proceedings of ICLR_.
  * Pei et al. (2020) Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. 2020.  Geom-GCN: Geometric graph convolutional networks. In _Proceedings of ICLR_.
  * Perozzi et al. (2014) Bryan Perozzi, Rami Al-Rfou, and Steven Skiena. 2014\.  Deepwalk: Online learning of social representations. In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_. 701–710.
  * Raghu et al. (2017) Maithra Raghu, Ben Poole, Jon Kleinberg, Surya Ganguli, and Jascha Sohl-Dickstein. 2017.  On the expressive power of deep neural networks. In _international conference on machine learning_. 2847–2854.
  * Rahman et al. (2020) Md Khaledur Rahman, Majedul Haque Sujon, and Ariful Azad. 2020\.  Force2Vec: Parallel force-directed graph embedding. In _2020 IEEE International Conference on Data Mining (ICDM)_. IEEE, 442–451.
  * Ribeiro et al. (2017) Leonardo FR Ribeiro, Pedro HP Saverese, and Daniel R Figueiredo. 2017.  struc2vec: Learning node representations from structural identity. In _Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining_. 385–394.
  * Rosenberg and Hirschberg (2007) Andrew Rosenberg and Julia Hirschberg. 2007.  V-Measure: A conditional entropy-based external cluster evaluation measure. In _Proceedings of EMNLP_.
  * Rozemberczki et al. (2021) Benedek Rozemberczki, Carl Allen, and Rik Sarkar. 2021.  Multi-Scale Attributed Node Embedding.  _Journal of Complex Networks_ 9, 2 (2021).
  * Saramäki et al. (2007) Jari Saramäki, Mikko Kivelä, Jukka-Pekka Onnela, Kimmo Kaski, and Janos Kertesz. 2007\.  Generalizations of the clustering coefficient to weighted complex networks.  _Physical Review E_ 75, 2 (2007), 027105.
  * Suresh et al. (2021) Susheel Suresh, Vinith Budde, Jennifer Neville, Pan Li, and Jianzhu Ma. 2021.  Breaking the Limit of Graph Neural Networks by Improving the Assortativity of Graphs with Local Mixing Patterns.  _Proceedings of KDD_ (2021).
  * Tang et al. (2015) Jian Tang, Meng Qu, Mingzhe Wang, Ming Zhang, Jun Yan, and Qiaozhu Mei. 2015\.  LINE: Large-scale information network embedding. In _Proceedings of WWW_.
  * Van der Maaten and Hinton (2008) Laurens Van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-SNE.  _Journal of machine learning research_ 9, 11 (2008).
  * van Dongen (2000) S. van Dongen. 2000\.  _Graph clustering by flow simulation_.  Ph. D. Dissertation. Utrecht University.
  * Van Dongen (2000) Stijn Van Dongen. 2000\.  Graph clustering by flow simulation.  _PhD thesis, University of Utrecht_ (2000).
  * Veličković et al. (2018) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. 2018.  Graph attention networks. In _Proceedings of ICLR_.
  * Wu et al. (2019) Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger. 2019\.  Simplifying graph convolutional networks. In _International conference on machine learning_. PMLR, 6861–6871.
  * Wu et al. (2020) Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip. 2020\.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ (2020).
  * Yang et al. (2020) Yiding Yang, Zunlei Feng, Mingli Song, and Xinchao Wang. 2020\.  Factorizable graph convolutional networks.  _Proceedings of NeurIPS_ (2020).
  * Zeng et al. (2020) Hanqing Zeng, Hongkuan Zhou, Ajitesh Srivastava, Rajgopal Kannan, and Viktor Prasanna. 2020\.  GraphSAINT: Graph sampling based inductive learning method. In _Proceedings of ICLR_.
  * Zhu et al. (2020) Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra. 2020\.  Beyond homophily in graph neural networks: Current limitations and effective designs.  _Advances in Neural Information Processing Systems_ (2020).

##  Appendix A Experimental Details

Software Versions. We use Python programming language v3.7.3 to implement and
perform all of our experiments. We use several python packages for our
implementation which are outlined as follows: PyTorch: 1.7.0, PyTorch
Geometric (Fey and Lenssen, 2019) version 1.6.1 along with PyTorch Sparse
0.6.8, PyTorch Scatter 2.0.5, PyTorch Cluster 1.5.8, PyTorch Spline Conv
1.2.0, NetworkX: 2.2, scikit-learn: 0.23.2, Matplotlib: 3.0.3, python-louvain:
0.13

Implementation Details. The proposed MarkovGCN model has been implemented in
PyG framework111https://github.com/rusty1s/pytorch_geometric. We employ both
$torch.mm$ and $torch\\_sparse.spspmm$ functions to perform dense and sparse
matrix multiplication operation, respectively, on line 4 of Algorithm 1. We
observe very little advantage of using the sparse matrix. Once MCL steps are
done, we select $L$ number of matrices from $k$ matrices based on the type of
MarkovGNN variant. Each convolutional layer of MarkovGNN has a different graph
structure which is different from GCN. Otherwise, similar types of operations
are performed in both forward and backward propagation as GCN using $GCNConv$
class. In $l$-th convolutional layer, we pass all edges of $A^{(l)}$ in
MarkovGCN whereas GCN uses edges of the original matrix $A$. Similar to other
GNN models, we employ the negative log-likelihood loss function and Adam
optimizer of PyTorch package. Instead of the threshold-based pruning strategy
in Algorithm 1, we have also implemented a sparse $k$-selection technique
where we keep top-$k$ weighted entries in each row of the stochastic matrix
after each iteration of MCL. However, we have found from empirical results
that this method is not very effective and thus we do not report here. To
conduct the experiment for layer-wise contribution of Fig. 7 (b), we set the
$\mathbf{weight.requires\\_grad}$ and $\mathbf{bias.requires\\_grad}$
parameters of the corresponding layer to _True_ while setting other model
parameters to _False_.

More Details on Dataset. We have collected USAir network from struc2vec
(Ribeiro et al., 2017) which is publicly
available222https://github.com/leoribeiro/struc2vec. Email network is
available at Stanford SNAP dataset333https://snap.stanford.edu/data/email-Eu-
core.html. Chameleon and Squirrel networks (Rozemberczki et al., 2021) have
been collected from FAGCN (Bo et al., 2021) which are publicly
available444https://github.com/bdy9527/FAGCN. PPI network has been collected
from SuiteSparse website555https://sparse.tamu.edu/Pajek/yeast. Cora and
Citeseer networks are available in PyG datasets666https://pytorch-
geometric.readthedocs.io/en/latest/modules/datasets.html. The last accessed
date for all datasets is Oct. 20, 2021.

Table 5. Parameters used to run other methods. We set embedding dimension to
64 for all methods. We run all GNNs in PyTorch-Geometric framework using Adam
Optimizer. Methods:Parameters
---
GCN (Kipf and Welling, 2017): number of convolutional layers = 2, ReLU
activation function, learning rate = 0.01
GAT (Veličković et al., 2018): number of convolutional layers = 2, head
features = 8, ELU activation function, dropout rate = 0.6, learning rate =
0.005
Cluster-GCN (Chiang et al., 2019): number of convolutional layers = 2, number
of partition = 32 and 64, ReLU activation function, dropout rate = 0.5,
learning rate = 0.005
GDC (Klicpera et al., 2019): column-wise normalization, PPR diffusion with
alpha=0.05, Top-K sparsification with k = 128, number of convolutional layers
= 2, ReLU activation function, learning rate = 0.01
GraphSAINT (Zeng et al., 2020): edge sampling approach, batch size = 128,
number of steps = 5, sample coverage = 15, number of convolutional layers = 3,
linear layer = 1, ReLU activation function, dropout rate = 0.2, learning rate
= 0.01
Table 6. Parameters used to run different variants of MarkovGCN in PyTorch-Geometric framework (Fig. 3, Fig. 4, and Table 3). We set embedding dimension to 64 for all methods. For all variants, we use Negative Log-Likelihood Loss Function and Adam Optimizer. $lrate$ \- Learning rate, $droprate$ \- dropout rate, $nlayers$ \- the number of convolutional layers. Datasets | MarkovGCN Variants | Parameters for Markov Clustering Steps | Parameters for Convolutional Layers
---|---|---|---
\- -_markov_agg_ | Threshold ($\theta$) | Inflation ($i$) | Row-Stochastic. | LeakyRelu | $lrate$ | $droprate$ | $nlayers$
USAir | False | 0.1 | 1.6 | FALSE | FALSE | 0.01 | 0.5 | 4
True | 0.09 | 1.5 | FALSE | FALSE | 0.01 | 0.5 | 4
Email | False | 0.3 | 1.5 | FALSE | FALSE | 0.01 | 0.4 | 3
True | 0.26 | 1.5 | FALSE | FALSE | 0.01 | 0.3 | 3
PPI | False | 0.75 | 1.7 | TRUE | FALSE | 0.01 | 0.1 | 3
True | 0.35 | 1.8 | TRUE | FALSE | 0.05 | 0.5 | 2
Chameleon | False | 0.06 | 1.8 | FALSE | FALSE | 0.01 | 0.7 | 2
True | 0.2 | 1.5 | TRUE | FALSE | 0.01 | 0.5 | 3
Squirrel | False | 0.2 | 1.5 | FALSE | FALSE | 0.005 | 0.2 | 2
True | 0.05 | 1.5 | FALSE | FALSE | 0.01 | 0.25 | 6
Cora | False | 0.28 | 1.5 | FALSE | FALSE | 0.01 | 0.9 | 2
True | 0.03 | 1.3 | FALSE | FALSE | 0.01 | 0.8 | 2
Citeseer | False | 0.11 | 1.6 | FALSE | FALSE | 0.01 | 0.5 | 2
True | 0.11 | 1.5 | FALSE | FALSE | 0.01 | 0.6 | 2
Actor | False | 0.35 | 1.5 | FALSE | FALSE | 0.01 | 0.5 | 3
True | 0.3 | 1.5 | FALSE | FALSE | 0.01 | 0.5 | 3
Table 7. The contingency table for the Chameleon network. G.T. - Ground Truth, $T_{i}$’s are the number of true labels and $P_{i}$’s represent the number of predicted labels. $(i,j)$th entry denotes the number of vertices in common between $T_{i}$ and $P_{j}$. The diagonal shaded cells represent the true positives for corresponding classes.  G.T. | GCN | GDC | MarkovGCN
---|---|---|---
| P1 | P2 | P3 | P1 | P2 | P3 | P1 | P2 | P3
T1 | 770 | 110 | 41 | 780 | 105 | 36 | 837 | 55 | 29
T2 | 75 | 527 | 98 | 86 | 543 | 71 | 94 | 551 | 55
T3 | 56 | 81 | 519 | 61 | 115 | 480 | 24 | 78 | 554

Classification, Clustering and Visualization using GNNs. We train all GNN
models for 200 epochs. After each epoch, we compute the validation and test
accuracy with the so far trained model. We choose the best value of the test
classification accuracy for all GNN models. The graph clustering task requires
labels for all vertices to avoid any bias in the experiment. Thus, we infer
the labels for all vertices using the trained model which are then compared
with the ground truth labels to compute performance metrics. For the
visualization task, we infer the 64-dimensional embedding from the last layer
of the GNNs, then apply t-SNE to generate 2D layouts using the default
parameters and plot them using $matplotlib$ package.

Hyper-parameters of MarkovGCN. We analyze the hyper-parameters of MarkovGCN
extensively using a grid search technique. For other MarkovGNNs, we only use
different graphs in different convolutional layers where other parameters
remain the same. In MarkovGCN, we analyze three parameters in the Markov
process and four parameters in convolutional layers. For Markov process, we
use the following strategies:

  * •

Set threshold ($\theta$) in the ranges $\\{0.001\ldots 0.009\\}$,
$\\{0.01\ldots 0.09\\}$, and $\\{0.1\ldots 0.9\\}$, and then increase by step
size $0.0005$, $0.005$, and $0.05$, respectively.

  * •

Set inflation ($r$) in the range $\\{1.1,\ldots 3.3\\}$ and increase by step
$0.2$. We do not observe much variation in performance when $r>2$.

  * •

Consider either row-stochastic or column-stochastic matrices in the Markov
process. By default, MarkovGCN uses column-stochastic matrices. If row-
stochastic parameter is set to _TRUE_ , then it uses row-stochastic matrices.

We explore the parameters for convolutional layers by setting activation
function as either LeakyReLU or ReLU, learning rate of Adam optimizer in the
ranges $\\{0.01\ldots 0.08\\}$ and $\\{0.001\ldots 0.009\\}$ increasing the
step size by $0.01$ and $0.001$, respectively. We set dropout rate in the
range $\\{0.1\ldots 0.9\\}$ increasing the step size by $0.1$ and the number
of convolutional layers from 2 to 10. Section 5.6 discusses the impact for
some of these parameters on MarkovGCN’s performance.

Figure 8. Sensitivity of $\alpha$ parameter for a homophilic and a
heterophilic networks.

Sensitivity of $\alpha$. In Fig. 8, we show the sensitivity of $\alpha$
parameter of Eqn. 4 using the Citeseer and Chameleon networks. We observe that
the homophilic Citeseer network tends to have a lower value of $\alpha$
whereas, the heterophilic Chameleon network tends to have a higher value to
produce better accuracy.

##  Appendix B Hyper Parameters

We use PyG implementations for other methods. The I/O part is the same for all
GNN methods. To run a model, we need to set some parameters. We collect those
values of parameters from the papers. If no information is provided, then we
use the default parameters. In Table 5, we report the parameters for other
methods.

We analyze the parameter sensitivity for different variants of MarkovGCN. We
use a grid search approach to search optimal value as described in Sec. A.
After an extensive set of experiments, we pick a suitable set of parameters
for the MarkovGCN models. We report the values of different parameters in
Table 6 that produced better results in our experiments. These sets of
parameters have been used to generate the results in Fig. 3, Fig. 4, and Table
3.

Table 8. Adjusted Rand Index (ARI) and V-Measure for DeepWalk and struc2vec for different graphs. Graphs | ARI | V-Measure
---|---|---
| DeepWalk | struc2vec | DeepWalk | struc2vec
USAir | 0.094 | 0.179 | 0.150 | 0.243
Email | 0.472 | 0.010 | 0.704 | 0.237
PPI | 0.030 | 0.033 | 0.107 | 0.068
Chameleon | 0.098 | 0.035 | 0.155 | 0.106
Squirrel | 0.013 | 0.019 | 0.031 | 0.025
Cora | 0.352 | 0.022 | 0.457 | 0.101
Citeseer | 0.205 | 0.011 | 0.212 | 0.066

##  Appendix C Additional results

Table 7 shows the contingency table for the Chameleon network where the
$(i,j)$th entry denotes the number of vertices in common between the ground
truth class $T_{i}$ and the predicted class $P_{j}$. Diagonal entries in Table
7 denote true positive pairs of vertices stay together in both ground truth
and predicted clusterings. These results are consistent with our hypothesis
that MarkovGNN predicts the clustering patterns if available in the input
graph.


# Fast interactive web-based data visualizer of panoramic spectroscopic
surveys

Ivan Katkov,1,2 Damir Gasymov,2,3 Joseph D. Gelfand,1 Viktoria Toptun,2,3
Kirill Grishin,2 Igor Chilingarian4,2 Anastasia Kasparova,2 Vladislav
Klochkov2,3 Evgenii Rubtsov,2,3 Vladimir Goradzhanov2,3

###### Abstract

Panoramic IFU spectroscopy is a core tool of modern observational astronomy
and is especially important for galaxy physics. Many massive IFU surveys, such
as SDSS MaNGA (10k targets), SAMI (3k targets), Califa (600 objects), Atlas3D
(260 objects) have recently been released and made publicly available to the
broad astronomical community. The complexity and massiveness of the derived
data products from spectral cubes makes visualization of the entire dataset
challenging, but nevertheless very important and crucial for scientific
output.

Based on our past experience with visualization of spectral and imaging data
built in the frame of the VOxAstro Initiative projects, we are now developing
online web service for interactive visualizing spectroscopic IFU datasets
(ifu.voxastro.org). Our service will provide a convenient access and
visualization tool for spectral cubes from publicly available surveys (MaNGA,
SAMI, Califa, Atlas3D) and results of their modeling, as well as maps of
parameters derived from cubes, implementing the connected views concept. Here
we describe the core components and functionality of the service, including
REST API implementation on top of the Django+Postgres backend as well as a
fast and responsive user interface built using the modern Vue.js-based
framework Quasar.

1New York University Abu Dhabi, Abu Dhabi, United Arab Emirates;
<EMAIL_ADDRESS>
2Sternberg Astronomical Institute, M.V. Lomonosov Moscow State University,
Moscow, Russia

3Department of Physics, M.V. Lomonosov Moscow State University, Moscow, Russia

4Smithsonian Astrophysical Observatory, Cambridge, MA, USA

##  1 Introduction

Massive spectroscopic surveys of galaxies provide extremely important
information for studying galaxy evolution such as stellar and ionized gas
kinematics, star formation history imprinted in the stellar population
properties. Integral Field Spectroscopy (IFS) is a powerful observational
technique that allows spatially resolved measurements of galaxy properties.

Here we present an online web service (hereafter IFU Visualiser) for
interactive visualization of spectroscopic IFU datasets
(https://ifu.voxastro.org), providing a convenient access and visualization
tool for spectral cubes from publicly available surveys (MaNGA, SAMI, Califa,
Atlas3D). The service is inspired by our first prototype
http://manga.voxastro.org and the official Marvin web application released by
the MaNGA Team (Cherinka et al. 2019).

##  2 Data sets and service functionality

To date (Nov 2021) IFU Visualiser involves data arrived from SDSS DR16 MaNGA
Survey (Bundy et al. 2015), SAMI DR3 (Croom et al. 2012), Califa DR3 (Sánchez
et al. 2012), Atlas3D (Cappellari et al. 2011). Table 1 shows the size of the
datasets available in the IFU Visualiser service and the cross-matching
between surveys.

Table 1.: IFU Visualiser data sets | MaNGA DR17 | MaNGA DR16 | SAMI | Califa | Atlas3D
---|---|---|---|---|---
MaNGA DR17 | 11236 | 4824 | 106 | 40 | 0
MaNGA DR16 | 4824 | 4824 | 15 | 20 | 0
SAMI | 106 | 15 | 3426 | 1 | 1
Califa | 40 | 20 | 1 | 667 | 23
Atlas3D | 0 | 0 | 1 | 23 | 260

IFU Visualiser largely inherits the general architecture and functionality of
another VOxAstro project RCSEDv2 (Chilingarian et al. 2017; Toptun 2021;
Klochkov 2021). Briefly, IFU Visualiser provides two major ingredients: search
tool for spectral cubes and interactive visualiser of individual cubes and
two-dimensional maps of parameters extracted from spectra. Query language
syntax described on the documentation page https://ifu.voxastro.org/docs where
several search query examples can be found. Fig. 1 shows a screenshot of the
individual Cube page.

Figure 1.:  Example of an individual page of the MaNGA spectral cube (Cube ID:
4295) showing the low-mass post-starburst galaxy GMP 4188 in the Coma cluster,
which will passively evolve into an ultra-diffuse galaxy (Grishin et al.
2021).

##  3 Implementation

All code, including the backend, frontend, and scripts to assembly the input
tables for the database, is available on the GitHub
repo.111https://github.com/voxastro/ifu-visualiser

###  3.1 Backend

PostgreSQL is used to manage access to the main and linked tables. The main
spectral cube table contains basic information such as cube identifier,
coordinates, field-of-view coverage, and the necessary survey identifiers
needed to connect related tables to any given cube. The linked tables are
taken from the corresponding survey and contain very diverse information about
cubes as well as galaxies. We chose PostgreSQL as the backend database because
of i) support for geometry queries thanks to q3c (Koposov & Bartunov 2006) and
pgSphere (Chilingarian et al. 2004) extensions; and ii) natural compatibility
with VO services (TAP, SSAP access) using GAVO DaCHS (Demleitner 2018), which
we plan to add soon to IFU Visualiser. The well-known Python-based Django REST
Framework (DRF)222https://www.django-rest-framework.org/ is used on top of the
PostgreSQL database to provide API endpoints. In addition to Django high
popularity, wide community, and excellent documentation, a decisive reason for
using this framework is the large number of batteries. One of them is `django-
custom-query`,333https://github.com/ivan-katkov/django-custom-query which is
designed to create simplified, but still very powerful SQL-like queries (e.g.
`survey=sami OR atlas_name~4551`). We also modified the original library and
incorporated coordinates Cone-search queries using the q3c function
`q3c_radial_query()` (e.g. `cone(195, 28, 0.5)`). We also considered using the
GraphQL444https://graphql.org/ endpoint to provide data access, but after
simple performance tests with Graphene-Django library555https://docs.graphene-
python.org/projects/django/en/latest/ we found that DRF provides several times
faster response time overall. Another problem with using GraphQL lies on the
frontend side: the seeming simplicity of using, for example, with Apollo
client666https://apollo.vuejs.org/ preserves potential problems in managing
global application state. The Apollo client can be used as a state manager,
but we use the standard Vuex state manager for Vue.js applications. Using both
options leads to the "two sources of truth" problem. Switching completely to
Apollo requires writing awkward chunks of code and abandoning browser testing
tools not yet available for the Apollo client.

###  3.2 Frontend

We used the Vue.js-based frontend framework Quasar777https://quasar.dev/,
which provides an awesome infrastructure and a large collection of Vue
components built according to Material Design guidelines. Vue.js and the
Quasar framework have a very low threshold of entry and allow one to create
highly interactive user interfaces with very limited Javascript skills. In
addition to the Quasar component library, we used Aladin Lite for imaging and
Plotly.js888https://plotly.com/javascript/ library to create an interactive
spectral plots and maps of parameters.

##  4 Next development steps

To date (November 2021), only a small part of the planned functionality of the
IFU Visualiser has been developed. The next major steps for further
development will be:

  1. 1.

Regarding user interface development we will improve the spectral Cube page
and implementing a more flexible interaction with spectra and parameter maps.
Functionality will be similar to that we built in the frame of our first
prototype of such a service https://manga.voxastro.org/ made with another
technologies and uses only MaNGA data.

  2. 2.

We also plan to add uniform basic galaxy parameters such as masses,
luminosities, effective radii and colors for all galaxies available in the IFU
Visualiser.

  3. 3.

Another important milestone is to analyze the entire collection of spectral
cubes to ensure a homogeneous set of data products extracted from the cubes
using full spectral fitting technique NBursts, similar to what we do in the
RCSED2 project but for IFU data.

Acknowledgments. This project is supported by the Russian Science Foundation
grant 21-72-00036. Authors also thank the Interdisciplinary Scientific and
Educational School of Moscow University “Fundamental and Applied Space
Research”. DG, KG, VT, EG are also thankful to the ADASS XXXI organizing
committee for providing financial aid.

## References

  * Bundy et al. (2015) Bundy, K., et al. 2015, ApJ, 798, 7. 1412.1482
  * Cappellari et al. (2011) Cappellari, M., et al. 2011, MNRAS, 413, 813. 1012.1551
  * Cherinka et al. (2019) Cherinka, B., et al. 2019, AJ, 158, 74. 1812.03833
  * Chilingarian et al. (2004) Chilingarian, I., et al. 2004, in Astronomical Data Analysis Software and Systems (ADASS) XIII, edited by F. Ochsenbein, M. G. Allen, & D. Egret, vol. 314 of Astronomical Society of the Pacific Conference Series, 225
  * Chilingarian et al. (2017) Chilingarian, I. V., et al. 2017, ApJS, 228, 14. 1612.02047
  * Croom et al. (2012) Croom, S. M., et al. 2012, MNRAS, 421, 872. 1112.3367
  * Demleitner (2018) Demleitner, M. 2018, DaCHS: Data Center Helper Suite. 1804.005
  * Grishin et al. (2021) Grishin, K. A., et al. 2021, Nature Astronomy. 2111.01140
  * Klochkov (2021) Klochkov, V. 2021, in ADASS XXX, edited by J.-E. Ruiz, & F. Pierfederici (San Francisco: ASP), vol. TBD of ASP Conf. Ser., 999 TBD
  * Koposov & Bartunov (2006) Koposov, S., & Bartunov, O. 2006, in Astronomical Data Analysis Software and Systems XV, edited by C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, vol. 351 of Astronomical Society of the Pacific Conference Series, 735
  * Sánchez et al. (2012) Sánchez, S. F., et al. 2012, A&A, 538, A8. 1111.0962
  * Toptun (2021) Toptun, V. 2021, in ADASS XXX, edited by J.-E. Ruiz, & F. Pierfederici (San Francisco: ASP), vol. TBD of ASP Conf. Ser., 999 TBD