# Community detection in multiplex networks based on orthogonal nonnegative
matrix tri-factorization

Meiby Ortiz-Bouza, Selin Aviyente  M. Ortiz-Bouza and S. Aviyente were with
the Department of Electrical and Computer Engineering, Michigan State
University, East Lansing, MI, 48824.
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Networks are commonly used to model complex systems. The different entities in
the system are represented by nodes of the network and their interactions by
edges. In most real life systems, the different entities may interact in
different ways necessitating the use of multiplex networks where multiple
links are used to model the interactions. One of the major tools for inferring
network topology is community detection. Although there are numerous works on
community detection in single-layer networks, existing community detection
methods for multiplex networks mostly learn a common community structure
across layers and do not take the heterogeneity across layers into account. In
this paper, we introduce a new multiplex community detection method that
identifies communities that are common across layers as well as those that are
unique to each layer. The proposed method, Multiplex Orthogonal Nonnegative
Matrix Tri-Factorization, represents the adjacency matrix of each layer as the
sum of two low-rank matrix factorizations corresponding to the common and
private communities, respectively. Unlike most of the existing methods which
require the number of communities to be pre-determined, the proposed method
also introduces a two stage method to determine the number of common and
private communities. The proposed algorithm is evaluated on synthetic and real
multiplex networks, as well as for multiview clustering applications, and
compared to state-of-the-art techniques.

###### Index Terms:

Multiplex Networks, Community Detection, Nonnegative Matrix Tri-factorization,
Eigengap, Low-Rank Structure

##  I Introduction

Complex networks are usually used to represent many real world systems,
ranging from social to biological ones [1], where the different agents and
their relations are represented as the nodes and edges of the network,
respectively. Traditional network models employ simple graphs, where there is
a single edge between any two nodes. Thus, these models cannot capture
multiple modes of interaction that may exist between the nodes. Recently,
multiplex networks that represent multiple modes of interaction have been
proposed. A multiplex network is a multilayer network where all layers share
the same set of nodes with different topologies [2]. Multiplex networks have
been used to model a variety of complex systems, including living organisms,
human societies, transportation systems and critical infrastructures [3, 4].

Community detection is an important tool in network analysis, where
communities are defined as groups of nodes that are more densely connected to
each other than they are to the rest of the network. Most of the existing work
on community detection [5] focuses on single layer networks. Community
detection methods for multiplex networks [6] can be grouped into three main
classes. The first class of methods merges the layers in a multiplex network,
using a flattening algorithm, then apply single-layer community detection to
the aggregated network [7, 8, 9]. While these methods are computationally
efficient, they can only identify communities that are common across all
layers. Moreover, due to the flattening process, some spurious communities may
emerge. The second class of methods applies community detection to each layer
individually and then merges the results [10, 11, 12]. These methods include
nodes in the same community only when they are part of the same community in
at least one layer. Finally, the third class of methods operates directly on
the multiplex network model [13, 14, 15, 16, 17, 18].

Existing multiplex community detection approaches typically assume that the
community structure is the same across layers and find the partition that best
fits all layers. Thus, they do not differentiate between communities that are
common across layers from those that are unique to each layer. This is
particularly important for real world applications where the networks are
heterogeneous, and the different layers correspond to different modes of
interaction. For example, in social networks, a group of individuals may be
well connected via friendships on Facebook; however, this group of individuals
will likely not work at the same company. Thus, in a situation like this, a
given community will only be present in a subset of the layers, and different
communities may be present in different subsets of layers.

In this work, common communities are defined as communities that are observed
in more than one layer, i.e., communities that are common across any subset of
two or more layers, and private communities as communities that are unique to
each layer. The problem of detecting common and private communities is then
formulated using a novel framework titled Multiplex Orthogonal Nonnegative
Matrix Trifactorization (MX-ONMTF). In the proposed framework, each layer’s
adjacency matrix is represented as the sum of two low-rank matrix
factorizations corresponding to the common and private communities,
respectively. The resulting optimization problem is solved using an iterative
multiplicative update algorithm. The proposed approach also addresses the
problem of determining the number of communities. Unlike most existing work,
where the number of communities is determined through a greedy search, in this
paper, a two-step approach is proposed. The proposed algorithm is first
evaluated on synthetic benchmark multiplex networks with different numbers of
layers, nodes, communities, noise levels, and inter-layer dependency
probability. Next, the proposed method is applied to real networks including
social and biological networks. Finally, the algorithm is evaluated for
multiview clustering task, where the communities across all layers are assumed
to be the same.

This paper extends our prior work [19] where the community detection problem
is formulated for two-layer multiplex networks. This work introduces
significant new contributions. First, MX-ONMTF is generalized for multiplex
networks with $L$ layers where the common communities can be observed for a
subset of layers. Second, a new two stage method is introduced to determine
the number of common communities such that communities that are common across
two or more layers can be identified. A detailed theoretical analysis of the
convergence of the algorithm and recovery guarantees is also provided.
Finally, this paper includes an extensive evaluation of the method on
synthetic and real networks, and multiview clustering applications.

The rest of the paper is organized as follows. Section II presents a summary
of related works. Section III provides background on community detection,
multiplex networks, and orthogonal nonnegative matrix tri-factorization.
Sections IV and V present the proposed multiplex community detection algorithm
and its convergence analysis. Section VI establishes the theoretical
properties of the algorithm, while Section VII illustrates results on both
simulated and real networks. Finally, Section VIII provides conclusions and
discussion on future work.

##  II Related works

The method proposed in this paper belongs to the third class of algorithms,
which operate directly on the multiplex network model. There are different
types of algorithms that fall in this class: random walk, statistical
generative network models, label propagation, objective function optimization,
and Nonnegative Matrix Factorization (NMF).

Methods based on random walkers model the dynamic process on networks as
random walks where the process is more likely to persist on the vertices in
the same community and far less on the vertices in different communities. For
instance, LART [13] is initialized by assigning each node in each layer to its
own community. Hierarchical clustering is then used to merge nodes based on a
distance matrix. The partition with the highest multiplex modularity is
chosen. In [14], Infomap which is based on a compression of network flows is
proposed to identify communities within and across layers. However, Infomap
tends to assign each physical node across layers to the same community, not
differentiating the topological differences across layers.

Statistical methods such as [20] use Weighted Stochastic Block Model (WSBM) to
detect common and private communities in heterogeneous weighted networks.
Although this method addresses the heterogeneity of networks across layers,
the method is limited to detecting only common communities that are shared by
all layers, ignoring communities that may be shared by only a subset or
different subsets of layers. In [21], authors propose a generative model and
an expectation maximization algorithm for community detection and link
prediction in multilayer networks. Although the method allows for different
connectivity patterns in each layer, the interdependence between layers is
only taken into account for link prediction, while the layers are assumed to
share a common community structure.

The third class of methods, Label Propagation Algorithms (LPA), is based on
the intuition that a label can become dominant in a densely connected group of
nodes but will have trouble crossing a sparsely connected region. In [22], an
LPA-based method for community detection in multidimensional networks is
proposed to identify communities and the subset of layers in which each of
these communities is observed, simultaneously. However, this algorithm fails
to detect communities that are private to each layer and communities that may
be common among a small number of layers.

The fourth type of multiplex community detection methods is based on defining
an objective function and identifying the community structure that
maximizes/minimizes the objective function. For example, Generalized Louvain
(GenLouvain) [15] uses an extended definition of modularity and is one of the
fastest methods for community detection in multiplex networks. As GenLouvain
assigns each node-layer tuple to its own community, it cannot identify common
communities across layers. More recently, multiobjective genetic and
evolutionary algorithms such as MultiMOGA [16] and MOEA/D-TS [23] have been
used to jointly maximize the modularity of each layer and the similarity
between the community structures across layers. These methods find a shared
community structure across all layers, not differentiating communities that
may be unique to each layer. In [24], extension of normalized cut to multiplex
networks is proposed by constructing a block Laplacian matrix with each block
corresponding to a layer. This method relies on selecting a parameter $\beta$
that controls the consistency of the community structure across different
layers.

The last class of methods is based on NMF which, because of its
interpretability and good performance, has been broadly used for community
detection in single-layer, multiplex, multilayer, and dynamic networks [25,
26, 27, 28]. In [29], Semi-Supervised joint Nonnegative Matrix Factorization
(S2-jNMF) is proposed for detecting the common communities across layers in a
multiplex network. A greedy search of dense subgraphs is performed and these
subgraphs are used as a priori information to create new adjacency matrices
for each layer. In [30], a two-step approach is proposed, where first a
nonnegative low dimensional feature representation of each layer is found
using one of the four different NMF models. These community structures are
then used to obtain a consensus community structure. Authors in [31] use NMF
for detecting communities in multiplex social networks, where both unifying
and coupling approaches are proposed. The unifying approach finds a common
community structure by aggregating all layers, while the coupling approach
finds mostly consistent community structures. Most of the aforementioned NMF
based methods find a common structure across all layers or for a majority of
layers and do not consider cases where common communities may be present in
different subsets of layers. Moreover, they do not detect private communities.
These methods also require that the number of communities is provided a
priori.

##  III Background

###  III-A Community Detection

Community detection for a single-layer network is the partitioning of a node
set $V$ as $\mathcal{C}=\\{C_{1},...,C_{K}\\}$ where $K$ is the number of
communities. One of the most popular algorithms for partitioning graphs is the
minimum cut method and its variants such as ratio cut and normalized cut [32].
In these methods, the network is partitioned such that the number of edges
between different communities is minimized.

Given a single layer graph, $\mathcal{G}=\\{V,E,\mathbf{A}\\}$, where $V$, $E$
and $\mathbf{A}$ are the set of nodes, edges, and adjacency matrix of the
graph, respectively, the min-cut problem aims to find a partition by
minimizing the following objective function:

$\frac{1}{2}\sum_{p=1}^{K}\mathbf{Cut}(C_{p},\bar{C}_{p}),$ (1)

where $\bar{C}_{p}$ is the complement of $C_{p}$ and
$\mathbf{Cut}(C_{p},\bar{C}_{p})=\sum_{i\in C_{p},j\in\bar{C}_{p}}A_{ij}$ for
two disjoint sets $C_{p}$ and $\bar{C}_{p}$. The min-cut problem is NP-hard.
However, it has been shown [32, 25] that spectral clustering and nonnegative
matrix factorization provide solutions to relaxed versions of the min-cut
problem.

In particular, spectral clustering solves the min-cut problem by embedding
each node in a lower dimensional subspace spanned by the eigenvectors of the
normalized Laplacian matrix $\mathbf{L}$, where
$\mathbf{L}=\mathbf{D}^{-1/2}(\mathbf{D-A})\mathbf{D}^{-1/2}$, with
$\mathbf{D}$ being the diagonal matrix defined as $D_{ii}=\sum_{j}A_{ij}$.
k-means clustering is then applied to this low-dimensional subspace spanned by
the eigenvectors. Authors in [25] show that NMF is equivalent to Laplacian
based spectral clustering. Normalized Cut using the normalized adjacency
matrix, $\tilde{\mathbf{A}}=\mathbf{D}^{-1/2}\mathbf{AD}^{-1/2}$, is
equivalent to the nonnegative matrix factorization problem
$\underset{\mathbf{H}\geq
0}{\operatorname{argmin}}=||\tilde{\mathbf{A}}-\mathbf{HH}^{\top}||^{2}$.

###  III-B Multiplex Networks

Multiplex networks can be represented using a finite sequence of graphs
$\\{\mathcal{G}_{l}\\}$, where $l\in\\{1,2,\ldots,L\\}$,
$\mathcal{G}_{l}=(V_{l},E_{l},\mathbf{A}_{l})$ [33]. $V_{l}$ is the set of
nodes in layer $l$ and $\mathbf{A}_{l}\in\mathbb{R}^{n\times n}$ is the
adjacency matrix for layer $l$. In this paper, we use undirected (symmetric)
weighted and binary adjacency matrices. For a weighted adjacency matrix,
$A_{l_{ij}}\in[0,1]$, and for a binary adjacency matrix,
$A_{l_{ij}}\in\\{0,1\\}$.

###  III-C Orthogonal Nonnegative Matrix Tri-Factorization

Nonnegative Matrix Factorization decomposes a nonnegative matrix
$\mathbf{W}\in\mathbb{R}^{n\times m}$ into the product of two low-rank
nonnegative matrices $\mathbf{V}\in\mathbb{R}^{n\times k}$ and
$\mathbf{U}\in\mathbb{R}^{m\times k}$, such that
$\mathbf{W}\approx\mathbf{VU}^{\top}$ and $k\ll n,m$. $\mathbf{V}$ and
$\mathbf{U}$ are found by solving the optimization problem

$\centering\underset{\textbf{V},\textbf{U}>0}{\operatorname{argmin}}||\mathbf{W}-\mathbf{VU}^{\top}||_{F}^{2}.\@add@centering$
(2)

In NMF-based community detection, $\mathbf{V}$ and $\mathbf{U}$ are the
community feature matrix and the community indicator matrix, respectively.
Adding orthogonality constraints improves the performance of NMF as
orthogonality and nonnegativity force each row of $\mathbf{V}$ ($\mathbf{U}$)
to have only one nonzero element which implies that each node belongs only to
one community. Orthogonal NMF has been broadly used in community detection
where $\mathbf{W}$ is usually the adjacency matrix and $k$ is the number of
communities [34, 35].

Figure 1: Illustration of the proposed community detection algorithm.

Different extensions of NMF, such as Symmetric Nonnegative Matrix Tri-
factorization [36], where $\mathbf{W}$ is approximated by
$\mathbf{USU}^{\top}$, have been proposed for community detection.
$\mathbf{U}$ contains community membership information while the symmetric
matrix $\mathbf{S}\in\mathbb{R}^{k\times k}$ provides more degrees of freedom
to the approximation. In this paper, we will use Orthogonal Nonnegative Matrix
Tri-factorization (ONMTF) to formulate the community detection problem in
multiplex networks.

##  IV Proposed Method (MX-ONMTF)

The proposed method, MX-ONMTF, models each layer’s adjacency matrix as a sum
of low-rank representations of common and private communities using Orthogonal
Nonnegative Matrix Tri-Factorization (ONMTF). Fig. 1 illustrates the overview
of the proposed algorithm for a multiplex network with $L$ layers and two
common communities.

###  IV-A Problem Formulation

In this paper, we define common communities as communities that are observed
in more than one layer.

Definition 1: An ideal common community $C$ in a multiplex network
$\\{\mathcal{G}_{l}\\}$ is defined as a subgraph with the same set of nodes
for a subset of layers $\textbf{m}\subseteq\\{1,2,\ldots,L\\}$, where
$|\textbf{m}|>1$. Mathematically, $C$ can be defined as

$\displaystyle\begin{split}C=&\\{(V^{C}_{l},E^{C}_{l}):V^{C}_{l}\subseteq
V_{l},E^{C}_{l}\subseteq V^{C}_{l}\times V^{C}_{l},V_{l}^{C}=V_{k}^{C}\\\
&l,k\in\textbf{m},\textbf{m}\subseteq\\{1,2,\ldots,L\\},|\textbf{m}|>1\\}.\\\
\end{split}$

Definition 2: A private community $C$ in a multiplex network
$\\{\mathcal{G}_{l}\\}$ is defined as any community that is not common across
at least two layers.

For a multiplex network with $L$ layers and adjacency matrices,
$\mathbf{A}_{l}\in\mathbb{R}^{n\times n}$, $l\in\\{1,2,\ldots,L\\}$, we model
each layer’s adjacency matrix in terms of common and individual communities
using ONMTF. The resulting objective function can be formulated as

$\displaystyle\underset{\mathbf{H}\geq 0,\mathbf{H}_{l}\geq
0,\mathbf{S}_{l}\geq 0,\mathbf{G}_{l}\geq
0}{\operatorname{argmin}}\sum_{l=1}^{L}||\mathbf{A}_{l}-\mathbf{HS}_{l}\mathbf{H}^{\top}-\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}||^{2}_{F}$
(3) $\displaystyle\mbox{s.t\hskip
2.84544pt}\mathbf{H}^{\top}\mathbf{H}=\mathbf{I},\mathbf{H}_{l}^{\top}\mathbf{H}_{l}=\mathbf{I},\mbox{with
\hskip 2.84544pt}l\in\\{1,2,\ldots,L\\},$

where $\mathbf{H}\in\mathbb{R}^{n\times k_{c}}$ and
$\mathbf{H}_{l}\in\mathbb{R}^{n\times k_{p_{l}}},l\in\\{1,2,\ldots,L\\}$ are
the community membership matrices corresponding to the common and private
communities, respectively, and $\mathbf{S}_{l}$ and $\mathbf{G}_{l}$ are
symmetric matrices. In this work, it is assumed that the $L$ layers have a
total of $k_{c}$ common communities and $k_{p_{l}}$ private communities in
each layer $l$. The goal is to simultaneously identify communities that are
common across any subset of two or more layers and communities that are unique
to each layer. Therefore, $\mathbf{H}$ will contain information for all common
communities.

###  IV-B Optimization solution

ONMTF optimization problem in (3) can be solved using a multiplicative update
algorithm (MUA) [36]. Multiplicative update algorithms for solving NMF
problems were introduced in [37], while solving NMTF with orthogonal
constraints was first addressed by [36]. In this paper, we follow their
approach to derive the multiplicative update rules for each variable.

To find the update rules for $\mathbf{H}$, $\mathbf{H}_{l}$, $\mathbf{S}_{l}$,
and $\mathbf{G}_{l}$, the following Lagrangian function with Lagrange
multipliers $\mathbf{\Lambda}$ and $\mathbf{\Lambda}_{l}$ is minimized:

$\displaystyle\begin{split}&\mathcal{L}(\mathbf{H},\mathbf{H}_{l},\mathbf{S}_{l},\mathbf{G}_{l})=\sum_{l=1}^{L}||\mathbf{A}_{l}-\mathbf{HS}_{l}\mathbf{H}^{\top}-\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}||^{2}_{F}\\\
&+tr(\mathbf{\Lambda(H}^{\top}\mathbf{H-I}))+\sum_{l=1}^{L}tr(\mathbf{\Lambda}_{l}(\mathbf{H}_{l}^{\top}\mathbf{H}_{l}-\mathbf{I})).\end{split}$
(4)

For updating $\mathbf{H}$, we find $\mathbf{\nabla_{H}\mathcal{L}}$ as

$\displaystyle\begin{split}\mathbf{\nabla_{H}\mathcal{L}}=&\sum_{l=1}^{L}(\mathbf{4HS}_{l}^{\top}\mathbf{H}^{\top}\mathbf{HS}_{l}+4\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}-4\mathbf{A}_{l}\mathbf{HS}_{l})\\\
+&\mathbf{4H\Lambda}.\end{split}$ (5)

Applying the KKT conditions $\mathbf{\nabla_{H}\mathcal{L}}=0$ and
$\mathbf{\nabla_{\Lambda}\mathcal{L}}=0$, we obtain:
(i)
$\mathbf{\Lambda}=\sum_{l=1}^{L}(\mathbf{-S}_{l}^{\top}\mathbf{S}_{l}-\mathbf{H}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}+\mathbf{H}^{\top}\mathbf{A}_{l}\mathbf{HS}_{l})$.
(ii) $\mathbf{H}^{\top}\mathbf{H=I}$.
Substituting (i) and (ii) in Eq. (5), we get

$\displaystyle\begin{split}\mathbf{\nabla_{H}\mathcal{L}}=&\sum_{l=1}^{L}(\mathbf{4H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}-4\mathbf{A}_{l}\mathbf{HS}_{l}+4\mathbf{HH}^{\top}\mathbf{A}_{l}\mathbf{HS}_{l}\\\
\mathbf{-4HH}^{\top}&\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}).\end{split}$
(6)

As discussed in [38], if the gradient of an error function,
$\mathbf{\varepsilon}$, is of the form
$\mathbf{\nabla\varepsilon=\nabla\varepsilon^{+}-\nabla\varepsilon^{-}}$,
where $\mathbf{\nabla\varepsilon^{+}>0}$ and
$\mathbf{\nabla\varepsilon^{-}>0}$, then the multiplicative update for
parameter $\mathbf{\Theta}$ has the form
$\mathbf{\Theta=\Theta\odot\frac{\nabla\varepsilon^{-}}{\nabla\varepsilon^{+}}}$.
It can be easily seen that the multiplicative update preserves the
nonnegativity of $\mathbf{\Theta}$, while $\mathbf{\nabla\varepsilon}=0$ when
the convergence is achieved. Following this procedure, from the gradient of
the error function in Eq. (6), we derive the following multiplicative update
rule for $\mathbf{H}$

$\mathbf{H}\leftarrow\mathbf{H}\odot\frac{\sum_{l=1}^{L}(\mathbf{A}_{l}\mathbf{HS}_{l}+\mathbf{HH}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l})}{\sum_{l=1}^{L}(\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}+\mathbf{HH}^{\top}\mathbf{A}_{l}\mathbf{HS}_{l})},$
(7)

where the multiplication and division are performed element-wise and both
numerator and denominator are positive. Similarly, we obtain the following
update rules for $\mathbf{H}_{l}$, $\mathbf{S}_{l}$, and $\mathbf{G}_{l}$, for
each $l\in\\{1,2,\ldots,L\\}$:

$\mathbf{H}_{l}\leftarrow\mathbf{H}_{l}\odot\frac{\mathbf{A}_{l}\mathbf{H}_{l}\mathbf{G}_{l}+\mathbf{H}_{l}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}^{\top}\mathbf{H}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}}{\mathbf{HS}_{l}^{\top}\mathbf{H}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}^{\top}+\mathbf{H}_{l}\mathbf{H}_{l}^{\top}\mathbf{A}_{l}\mathbf{H}_{l}\mathbf{G}_{l}},\\\
$ (8)

$\mathbf{S}_{l}\leftarrow\mathbf{S}_{l}\odot\frac{\mathbf{H}^{\top}\mathbf{A}_{l}\mathbf{H}}{\mathbf{H}^{\top}\mathbf{HS}_{l}\mathbf{H}^{\top}\mathbf{H}+\mathbf{H}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}},$
(9)

$\mathbf{G}_{l}\leftarrow\mathbf{G}_{l}\odot\frac{\mathbf{H}_{l}^{\top}\mathbf{A}_{l}\mathbf{H}_{l}}{\mathbf{H}_{l}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}_{l}+\mathbf{H}_{l}^{\top}\mathbf{H}\mathbf{S}_{l}\mathbf{H}^{\top}\mathbf{H}_{l}}.$
(10)

Since NMF algorithms are initialized with random matrices, different runs
yield local minima. For this reason, we run the algorithm 50 times and report
the best results [39, 40]. As shown in Algorithm 1, for each random
initialization of $\mathbf{H}$, $\mathbf{H}_{l}$, $\mathbf{S}_{l}$, and
$\mathbf{G}_{l}$, the multiplicative update rules described in Eqs. (7)-(10)
are repeated for 1000 iterations or until convergence. We then select the
solution that yields the maximum value of the performance metric across the
different runs. For synthetic networks for which a ground truth is available,
Normalized Mutual Information (NMI) [41] is used. For networks without ground
truth, Modularity Density ($Q_{D}$) [42] is used as the performance metric.

Algorithm 1 MX-ONMTF

0: Adjacency matrices $\mathbf{A}_{l}$, $l\in\\{1,2,\ldots,L\\}$.

0: Community membership matrices $\mathbf{H}$, $\mathbf{H}_{l}$.

1: Use Algorithm 2 to find $k_{c}$ and $k_{p_{l}}$.

2: for r=1 to 50 do

3: Randomly initialize
$\mathbf{H},\mathbf{H}_{l},\mathbf{S}_{l},\mathbf{G}_{l}>0$

4: for 1000 iterations or until convergence do

5: update $\mathbf{H}$ according to Eq. (7)

6: update $\mathbf{H}_{l}$ for each $l\in\\{1,2,\ldots,L\\}$ according to Eq.
(8)

7: update $\mathbf{S}_{l}$ for each $l\in\\{1,2,\ldots,L\\}$ according to Eq.
(9)

8: update $\mathbf{G}_{l}$ for each $l\in\\{1,2,\ldots,L\\}$ according to Eq.
(10)

9: end for

10: for each layer $l$ do

11: Apply Algorithm 3 with $\mathbf{A}_{l}$, $\mathbf{H}$, and $k_{p_{l}}$ as
inputs to to find $\mathbf{H}_{c_{l}}$.

12: for each i do

13: $j^{*}\leftarrow argmax_{j}\mathbf{H}_{c_{l}}(i,j)$

14: if $\mathbf{H}_{c_{l}}(i,j^{*})>\mathbf{H}_{l}(i,j^{*})$ then

15: $idx(i)\leftarrow j^{*}$

16: else

17:
$idx(i)\leftarrow(argmax_{j}\mathbf{H}_{l}(i,j))+k_{c}+\sum_{n=1}^{l-1}k_{p_{n}}$

18: end if

19: end for

20: end for

21: Compute $\text{NMI}_{r}$ or $Q_{D_{r}}$.

22: end for

23: Choose the partition
$r^{*}=argmax_{r}\text{NMI}_{r}$($r^{*}=argmax_{r}Q_{D_{r}}$ ).

###  IV-C Number of communities

In most NMF-based community detection algorithms, the number of communities
($k$) is an input parameter. This problem is usually addressed by detecting
communities with different values of $k$ and selecting the one that gives the
solution with the best pre-determined performance metric, such as modularity
[43].

Algorithm 2 Finding $k_{l}$, $k_{c}$ and $k_{p_{l}}$, for
$l\in\\{1,2,\ldots,L\\}$ .

0: Adjacency matrices $\mathbf{A}_{l}$ for $l\in\\{1,2,\ldots,L\\}$.

0: Number of common communities $k_{c}$, total number of communities per layer
$k_{l}$, and number of private communities per layer $k_{p_{l}}$, for
$l\in\\{1,2,\ldots,L\\}$.

1: Let $\mathbf{L}^{null}_{l}$ be the normalized Laplacian the of Erdős–Rényi
null model.

2:
$\mathbf{L}^{null}_{l}=\mathbf{V}_{l}\mathbf{\Lambda}_{l}\mathbf{V}_{l}^{T}$

3: $\delta\leftarrow
quantile_{0.95}[max\\{|\lambda^{null}_{i}|-|\lambda^{null}_{i+1}|,i\geq 2\\}]$

4: $k_{l}\leftarrow min\\{k:|\lambda_{i}|-|\lambda_{i+1}|>\delta,\forall
i>k\\}$

5: Randomly initialize $\mathbf{U}_{l}\geq 0,\mathbf{S}_{l}\geq 0$.

6: for 1000 iterations or until convergence do

7: update $\mathbf{U}_{l}$ using
$\mathbf{U}_{l}=\mathbf{U}_{l}*\frac{(\mathbf{A}_{l}\mathbf{U}_{l}\mathbf{S}_{l})_{ij}}{(\mathbf{U}_{l}\mathbf{U}_{l}^{\top}\mathbf{A}_{l}\mathbf{U}_{l}\mathbf{S}_{l})_{ij}}$

8: update $\mathbf{S}_{l}$ using
$\mathbf{U}_{l}=\mathbf{U}_{l}*\frac{(\mathbf{U}_{l}^{\top}\mathbf{A}_{l}\mathbf{U}_{l})_{ij}}{(\mathbf{U}_{l}^{\top}\mathbf{U}_{l}\mathbf{S}_{l}\mathbf{U}_{l}^{\top}\mathbf{U}_{l})_{ij}}$

9: end for

10:
$\mathbf{X}=[\mathbf{U}_{1},\mathbf{U}_{2},\ldots,\mathbf{U}_{L}]^{\top}\in\mathbb{R}^{m\times
n}$, $m=\sum_{l=1}^{L}k_{l}$

11: $\mathbf{F}\leftarrow AgglomerativeHierarchicalClustering(\mathbf{X})$

12: $k_{c}=0$

13: for $i=2$ to $m-2$ do

14:
$d_{i}\leftarrow\frac{\mathbf{F}(i,3)-\mathbf{F}(i-1,3)}{\mathbf{F}(i-1,3)}$

15: if $d_{i}\geq 0.5$ then

16: if $max\\{\mathbf{F}(i,1),\mathbf{F}(i,2)\\}\leq m$ then

17: $k_{c}\leftarrow k_{c}+1$

18: else

19: $k_{c}\leftarrow k_{c}$

20: end if

21: else

22: $cut\leftarrow i$ and stop for

23: end if

24: end for

25:
$C\leftarrow$find$(1+\sum_{j=1}^{l-1}k_{j}\leq\mathbf{F}(1:cut,1:2)\leq\sum_{j=1}^{l}k_{j})$

26: $k_{p_{l}}\leftarrow k_{l}-|C|$, for $l\in\\{1,2,\ldots,L\\}$

In this paper, a two-step approach is proposed to determine the number of
communities per layer and the number of common communities. First, the number
of communities per layer ($k_{1}$, $k_{2}$,…, $k_{L}$), is found using the
eigengap rule [44]. Next, ONMTF is applied to each layer [36] and the low-rank
embedding matrices, $\mathbf{U}_{l}\in\mathbb{R}^{n\times k_{l}}$, are
obtained. Each element of the embedding matrices, $\mathbf{U}_{l}(i,j)$,
represents the likelihood of node $i$ belonging to community $j$. An
agglomerative hierarchical clustering algorithm using Euclidean distance is
applied on the rows of
$X=[\mathbf{U}_{1},\mathbf{U}_{2},\ldots,\mathbf{U}_{L}]^{\top}\in\mathbb{R}^{m\times
n}$, where $m=\sum_{l=1}^{L}k_{l}$, to obtain the number of common
communities. At each step of the algorithm, the two columns with the smallest
distance are aggregated, and the distances between the newly formed cluster
and the remaining ones are updated. A dendrogram like the one shown in Fig. 2
can be used to represent the different iterations of this algorithm. The $m$
leaves of the dendrogram correspond to the total number of communities across
the $L$ layers.

This agglomerative hierarchical clustering algorithm outputs a matrix
$\mathbf{F}\in\mathbb{R}^{m-1\times 3}$. The first two columns of
$\mathbf{F}(i,:)$ correspond to the labels of the two leaves of the dendrogram
that form cluster $m+i$ and the third column contains the distance between
these two leaves. This distance matrix $\mathbf{F}$ is used in Algorithm 2 to
determine the number of common communities, $k_{c}$, and the number of private
communities per layer, $k_{p_{l}}$. The algorithm iterates until the minimum
distance between any two clusters increases by more than $50\%$ of the minimum
distance from the previous iteration. Fig. 2 shows the dendrogram of the
hierarchical clustering of the columns of the embedding matrices
$\mathbf{U}_{1},\mathbf{U}_{2},\mathbf{U}_{3}$ of a 3-layer network with
$k_{1}=6$, $k_{2}=6$, and $k_{3}=5$ and the red line indicates where the
algorithm stops. For this example, $k_{c}=3$, $k_{p_{1}}=3$, $k_{p_{2}}=4$,
and $k_{p_{3}}=3$.

Figure 2: Example dendrogram illustrating the hierarchical clustering of the
columns of the embedding matrices
$\mathbf{U}_{1},\mathbf{U}_{2},\mathbf{U}_{3}$ for a 3-layer network with 3
common communities. The red line indicates where the algorithm stops.

###  IV-D Determining the common community labels for each layer

$\mathbf{H}\in\mathbb{R}^{n\times k_{c}}$ is the community membership matrix
corresponding to the common communities. Each row of $\mathbf{H}$ indicates
whether a particular node belongs to any of the $k_{c}$ common communities.
Algorithm 3 describes how the membership of the common communities is
determined across layers, i.e., which columns of $\mathbf{H}$ contain
information about the common communities present in layer $l$.

Algorithm 3 Identify the membership of $k_{c}$ common communities across
layers $\\{1,2,\ldots,L\\}$

0: Community membership matrices $\mathbf{H},\mathbf{H}_{l}$, adjacency
matrices $\mathbf{A}_{l}$, and number of private communities $k_{p_{l}}$, for
each layer $l\in\\{1,2,\ldots,L\\}$.

0: Layer specific community membership matrix $\mathbf{H}_{c_{l}}$ containing
information about the common communities present in layer $l$.

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

2: for each node $i$ do

3: $j^{*}\leftarrow argmax_{j}(\mathbf{[H,H}_{l}])_{ij}$

4: $idx(i)\leftarrow j^{*}$

5: end for

6: for $j=1$ to $k_{c}$ do

7: $\mathbf{H}(find(idx==j),j)=1$

8: $\mathbf{H}(find(idx\neq j),j)=0$

9: end for

10: for $j=1$ to $k_{c}$ do

11: $m=|\mathbf{H}(:,j)==1|$

12:
$\mathbf{B}_{0}\leftarrow\mathbf{A}_{l}\odot(\mathbf{H}(:,j)\mathbf{H}(:,j)^{\top})$

13: $\mathbf{B}_{1}\leftarrow\mathbf{A}_{l}\odot(\mathbf{1}_{n\times
n}-\mathbf{H}(:,j)\mathbf{H}(:,j)^{\top})$

14:
$q(j)\leftarrow\frac{\sum_{v,w}(\mathbf{B}_{1})_{vw}/(m(m-1))}{\sum_{v,w}(\mathbf{B}_{0})_{vw}/(m(n-m))}$

15: end for

16: $[q_{sorted},J]\leftarrow sort($q$,descending)$. $J$ contains the sorted
indices of the elements of $q$.

17: $pos\leftarrow J(1:k_{p_{l}})$.

18: $\mathbf{H}_{c_{l}}\leftarrow\mathbf{H}(:,[pos])$

19: end for

###  IV-E Time complexity

The time complexity of the proposed algorithm is mostly due to the
Multiplicative Updates Rules, Eqs. (7)-(10). The time complexity for the
product of two matrices, e.g., the product of a $m\times k$ matrix by a
$k\times n$ matrix, is $\mathcal{O}(mkn)$. Table II shows the time
complexities of Eqs. (7)-(10) and the total complexity, with
$l\in\\{1,2,\ldots,L\\}$, and $k=k_{c}+\sum_{l=1}^{L}k_{p_{l}}$.

$\mathbf{H}$ update (Eq. (7)) | $\mathcal{O}(n^{2}(k_{c}+\sum_{l=1}^{L}k_{p_{l}}))$
---|---
$\mathbf{H}_{l}$ update (Eq. (8)) | $\mathcal{O}(n^{2}(k_{c}+k_{p_{l}}))$
$\mathbf{S}_{l}$ update (Eq. (9)) | $\mathcal{O}(n^{2}k_{c})$
$\mathbf{G}_{l}$ update (Eq. (10)) | $\mathcal{O}(n^{2}k_{p_{l}})$
Total | $\mathcal{O}(n^{2}k)$
TABLE I: Computational complexity of updating each variable per iteration.

###  IV-F Storage complexity

The storage complexity of our algorithm is determined by the sizes of the
matrices $\mathbf{H}$, $\mathbf{H}_{l}$, $\mathbf{S}_{l}$, and
$\mathbf{G}_{l}$. It can be seen that the total storage complexity is
$\mathcal{O}(nk)$. For a multiplex network of size $n\times n\times L$, this
is a significant reduction in memory cost.

$\mathbf{H}$ | $\mathcal{O}(nk_{c})$
---|---
$\mathbf{H}_{l}$ | $\mathcal{O}(nk_{p_{l}})$
$\mathbf{S}_{l}$ | $\mathcal{O}({k_{c}}^{2})$
$\mathbf{G}_{l}$ | $\mathcal{O}({k_{p_{l}}}^{2})$
Total | $\mathcal{O}(nk)$
TABLE II: Storage complexity of each variable.

##  V Convergence Analysis

In this section, we will prove the convergence of the multiplicative update
rule defined by Eq. (7) using the auxiliary function approach. As the other
update rules are similar, we will not explicitly prove their convergence. We
first introduce the definition of auxiliary function as follows.

Definition 1: A function $Z(\mathbf{H},\mathbf{H}^{t})$ is called an auxiliary
function of $\mathbf{\mathcal{L}(H)}$ if it satisfies

$Z(\mathbf{H,H}^{t})\geq\mathbf{\mathcal{L}(H)}$ and
$Z(\mathbf{H,H})=\mathbf{\mathcal{L}(H)}$.

The auxiliary function is a useful concept because of the following lemma
which is proved in [37].

######  Lemma 1.

If $Z$ is an auxiliary function, then $\mathcal{L}$ is non-increasing under
the update

$\mathbf{H}^{t+1}=\underset{\mathbf{H}}{\operatorname{argmin}}Z(\mathbf{H,H}^{t})$.

######  Theorem 1.

Given $\mathbf{H}_{l}$, $\mathbf{S}_{l}$, and $\mathbf{G}_{l}$ the Lagrangian
function $\mathbf{\mathcal{L}(H)}$ is monotonically decreasing under the
update rule (7).

###### Proof.

For convenience, let $\mathcal{L}(h)$ denote the part of
$\mathcal{L}(\mathbf{H})$ dependent on $H_{ij}$. From Eq. (6) we have

$\displaystyle\begin{split}&\mathcal{L^{\prime}}(h)=\sum_{l=1}^{L}(\mathbf{4H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l}-4\mathbf{A}_{l}\mathbf{HS}_{l}\\\
&+4\mathbf{HH}^{\top}\mathbf{A}_{l}\mathbf{HS}_{l}\mathbf{-4HH}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l})_{ij}.\\\
\end{split}$

The second-order derivative of $\mathcal{L}(h)$ with respect to $h_{ij}$ is

$\displaystyle\begin{split}\mathcal{L^{\prime\prime}}(h)=\sum_{l=1}^{L}\\{4(\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top})_{ii}\mathbf{S}_{l_{jj}}-4\mathbf{A}_{l_{ii}}\mathbf{S}_{l_{jj}}\\\
+4[(\mathbf{H}^{\top}\mathbf{A}_{l}\mathbf{HS}_{l})_{ij}+h_{ij}(\mathbf{A}_{l}\mathbf{HS}_{l})_{ij}+(\mathbf{HH}^{\top}\mathbf{A}_{l})_{ii}\mathbf{S}_{l_{jj}}]\\\
-4[(\mathbf{H}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l})_{ij}+h_{ij}(\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{HS}_{l})_{ij}\\\
+(\mathbf{HH}^{\top}\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top})_{ii}\mathbf{S}_{l_{jj}}]\\}.\end{split}$

Let $h^{t}_{ij}$ denote the updated value of $h_{ij}$ after the $t$th
iteration, then the Taylor series expansion of $\mathcal{L}(h)$ at
$h^{t}_{ij}$ can be written as

$\displaystyle\begin{split}\mathcal{L}(h)=\mathcal{L}(h^{t}_{ij})+\mathcal{L^{\prime}}(h^{t}_{ij})(h-h^{t}_{ij})+\frac{1}{2}\mathcal{L^{\prime\prime}}(h_{ij}^{t})(h-h^{t}_{ij})^{2}.\end{split}$

Now, the key is to find an appropriate auxiliary function $Z(h,h_{ij}^{t})$.
We choose the following $Z(h,h_{ij}^{t})$ and prove in Appendix A, that it
satisfies the conditions to be an auxiliary function of $\mathcal{L}(h)$.

$\displaystyle\begin{split}&Z(h,h_{ij}^{t})=\mathcal{L}(h^{t}_{ij})+3\mathcal{L^{\prime}}(h^{t}_{ij})(h-h^{t}_{ij})+\\\
&\frac{3}{2}\frac{\sum\limits_{l=1}^{L}(4\mathbf{H}^{t}_{l}\mathbf{G}_{l}{\mathbf{H}_{l}^{t}}^{\top}\mathbf{H}^{t}\mathbf{S}_{l}+4\mathbf{H}^{t}{\mathbf{H}^{t}}^{\top}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}(h-h^{t}_{ij})^{2}.\end{split}$
(11)

According to Lemma 1, we must find the minimum of $Z(h,h_{ij}^{t})$ with
respect to $h$.

$\displaystyle\frac{\partial Z(h,h_{ij}^{t})}{\partial h}=$ $\displaystyle
3\mathcal{L^{\prime}}(h^{t}_{ij})+3\frac{\sum_{l=1}^{L}(4\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}(h-h^{t}_{ij})$
$\displaystyle+$ $\displaystyle
3\frac{\sum_{l=1}^{L}(4\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}(h-h^{t}_{ij})=0$

Replacing $\mathcal{L^{\prime}}(h^{t}_{ij})$ in the equation above and
canceling the common terms, we obtain

$\displaystyle\begin{split}\sum_{l=1}^{L}(-4\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l}\mathbf{-4H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}\\\
+\sum_{l=1}^{L}(4\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l}+4\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}\frac{h}{h_{ij}^{t}}=0.\end{split}$

Replacing $h$ by $h_{ij}^{t+1}$ we obtain the following update rule

$\displaystyle
h^{t+1}_{ij}=h^{t}_{ij}\frac{\sum_{l=1}^{L}(\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l}+\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{\sum_{l=1}^{L}(\mathbf{H}_{l}\mathbf{G}_{l}^{\top}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l}+\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}},$

which is the same as the update rule shown in Eq. (7). ∎

##  VI Recovery Guarantees

In this section, we will establish the theoretical properties of the proposed
community detection method. We want to determine if, by optimizing the
objective function, our algorithm will return good community structures. In
particular, we will investigate the consistency properties of the global
optimizer of the objective function under the multilayer stochastic blockmodel
(MLSBM). The optimization problem in (3) can be rewritten as

$\displaystyle\underset{\mathbf{H^{\prime}}_{l},\mathbf{F}_{l}}{\operatorname{argmin}}\sum_{l=1}^{L}||\mathbf{A}_{l}-\mathbf{H^{\prime}}_{l}\mathbf{F}_{l}\mathbf{H^{\prime}}_{l}^{\top}||^{2}_{F}\mbox{\hskip
5.69046pts.t\hskip
5.69046pt}\mathbf{H^{\prime}}_{l}^{\top}\mathbf{H^{\prime}}_{l}=\mathbf{I},$
(12)

where $\mathbf{F}_{l}$ is a block matrix defined as

$\mathbf{F}_{l}=\begin{bmatrix}\mathbf{S}_{l}&\bf{0}\\\
\bf{0}&\mathbf{G}_{l}\end{bmatrix}$

and $\mathbf{H^{\prime}}_{l}=[\mathbf{H}|\mathbf{H}_{l}]$ is the concatenation
of the community membership matrices of the common and private communities,
$\mathbf{H}\in\mathbb{R}^{n\times kc}$ and
$\mathbf{H}_{l}\in\mathbb{R}^{n\times k_{p_{l}}}$, respectively.

For the $n\times n\times L$ population adjacency tensor
$\mathcal{A}=\\{\mathbf{A}_{1},...,\mathbf{A}_{L}\\}$, we can define a
multiplex SBM $[\mathbf{Z},\Theta]$ as in [20], with each of the $L$ slices
$\mathbf{A}_{l}\in\mathbb{R}^{n\times n}$. The mulitplex SBM with parameters
$[\mathcal{Z}=\\{\mathbf{Z}_{1},...,\mathbf{Z}_{L}\\},\mathcal{B}=\\{\mathbf{\Theta}_{1},...,\mathbf{\Theta}_{L}\\}]$,
can be written in the matrix form as,

$\mathbb{E}(\mathbf{A}_{l})=\mathbf{Z}_{l}\mathbf{\Theta}_{l}\mathbf{Z}_{l}^{\top},$

with $\mathbf{\Theta}_{l}\in[0,1]^{(kc+k_{p_{l}})\times(kc+k_{p_{l}})}$ and
$\mathbf{Z}_{l}\in\\{0,1\\}^{n\times(kc+k_{p_{l}})}$ for each layer $l$.
$\mathbf{\Theta}_{l}$ is a block matrix defined as

$\mathbf{\Theta}_{l}=\begin{bmatrix}\mathbf{\Theta}_{c}&\bf{0}\\\
\bf{0}&\mathbf{\Theta}_{p_{l}}\end{bmatrix},$

with $\mathbf{\Theta}_{c}$ and $\mathbf{\Theta}_{p_{l}}$ being the affinity
probability matrices of the common and private communities, respectively.
$\mathbf{Z}_{l}=[\mathbf{Z}_{c}|\mathbf{Z}_{p_{l}}]$ is the concatenation of
the community membership matrices of the common and private communities,
$\mathbf{Z}_{c}\in\mathbb{R}^{n\times kc}$ and
$\mathbf{Z}_{p_{l}}\in\mathbb{R}^{n\times k_{p_{l}}}$, respectively.

To prove that our method can correctly recover the community assignments, we
propose the following lemma following the work in [45].

######  Lemma 2.

The optimization problem in (3) applied to $\mathcal{A}$ has
$\mathbf{H^{\prime}}_{l}=\mathbf{Z}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}$
and
$\mathbf{F}_{l}=(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}\mathbf{\Theta}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}$,
$l=1,...,L$, as the unique solution up to an orthogonal matrix, provided at
least one of the $\mathbf{\Theta}_{l}$ is full rank.

###### Proof.

To prove lemma 2, we can show that
$\mathbf{H^{\prime}}_{l}=\mathbf{Z}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}$,
$\mathbf{F}_{l}=(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}\mathbf{\Theta}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}$
is a solution to the optimization problem in (12). Substituting the solution
to (12), we have

$\centering\begin{split}\sum_{l=1}^{L}||\mathbf{A}_{l}-\mathbf{Z}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}&(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}\mathbf{\Theta}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}\mathbf{Z}_{l}^{\top}||^{2}_{F}\\\
&=\sum_{l=1}^{L}||\mathbf{A}_{l}-\mathbf{Z}_{l}\mathbf{\Theta}_{l}\mathbf{Z}_{l}^{\top}||^{2}_{F}\end{split}\@add@centering$

and, since
$\mathbb{E}(\mathbf{A}_{l})=\mathbf{Z}_{l}\mathbf{\Theta}_{l}\mathbf{Z}_{l}^{\top}$,
the value of this minimization objective function is 0, and
$\mathbf{H^{\prime}}_{l}^{\top}\mathbf{H^{\prime}}_{l}=(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}=(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{1/2}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}=\mathbf{I}$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Mean NMI over 100 realizations of (a)-(c) 3-layer, 4-layer and
5-layer benchmark networks, respectively, for the scenario with 2 common
communities across all layers; (d)-(f) 3-layer, 4-layer and 5-layer benchmark
networks, respectively, for the scenario with 3 common communities across
different subsets of layers. All networks are generated with 8 different
values of the mixing parameter $\mu$ and $n=256$.

Now, we need to show the uniqueness of this solution. By assumption, at least
one of the $\mathbf{\Theta}_{l}$ is full rank. For a non-singular matrix
$\mathbf{P}\in\mathbb{R}^{(k_{c}+k_{p_{l}})\times(k_{c}+k_{p_{l}})}$ we can
say that $\mathbf{H^{\prime}}_{l}\mathbf{P}$ and
$\mathbf{P}^{-1}\mathbf{\Theta}_{l}\mathbf{P}^{{-1}^{\top}}$ is also a
solution. Due to the orthogonality constraint, we must have
$(\mathbf{H^{\prime}}_{l}\mathbf{P})^{\top}\mathbf{H^{\prime}}_{l}\mathbf{P}=\mathbf{P}^{\top}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l}(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}\mathbf{P}=\mathbf{I}$,
which implies $\mathbf{P}^{\top}\mathbf{P}=\mathbf{I}$, and therefore the
solution is unique up to an orthogonal matrix. Moreover, since
$\mathbf{Q}^{-1/2}=(\mathbf{Z}_{l}^{\top}\mathbf{Z}_{l})^{-1/2}$ is a diagonal
matrix with positive elements and therefore invertible, we have that
$\mathbf{Z}_{i}\mathbf{Q}^{-1/2}=\mathbf{Z}_{j}\mathbf{Q}^{-1/2}$ implies
$\mathbf{Z}_{i}=\mathbf{Z}_{j}$. ∎

##  VII Experiments

###  VII-A Synthetic Multiplex Networks

####  VII-A1 Model description

Multiplex benchmark networks based on the model described in [46, 47] were
generated. The authors in [46] propose a two-step approach to generate
multilayer networks with a community structure. First, a multilayer partition
with the user-defined number of nodes in each layer, number of layers, and an
interlayer dependency tensor that specifies the desired dependency structure
between layers is generated. Next, for the given multilayer partition, edges
in each layer are generated following a degree-corrected block model [48]
parameterized by the distribution of expected degrees and a community mixing
parameter $\mu\in[0,1]$. The mixing parameter $\mu$ controls the modularity of
the network. When $\mu=0$, all edges lie within communities, whereas $\mu=1$
implies that edges are distributed independently. For multiplex networks, the
probabilities in the interlayer dependency tensor are the same for all pairs
of layers and are specified by $p\in[0,1]$. When $p=0$, the partitions are
independent across layers while $p=1$ indicates an identical partition across
layers.

In this paper, we extend the model described above to generate multiplex
benchmark networks with common and private communities. We first generate the
common communities by randomly selecting $n_{c}$ nodes across all layers and
setting the inter-layer dependency probability to $p_{1}$. For each common
community, we decide whether it exists in a particular layer or not. Next, we
independently generate the private communities for each layer with the
remaining nodes in that layer.We generated 100 different random realizations
of each multiplex network in order to report the average performance metric on
the experiments.

####  VII-A2 Evaluation

We compared the performance of our method to well-known multiplex community
detection algorithms. In particular, we compared with ONMTF applied to the
aggregated multiplex networks using the average of the adjacency matrices
(Aggregated Average), Spectral Clustering on Multi-Layer graphs (SC-ML) [12],
Generalized Louvain (GenLouvain) multilayer community detection algorithm [49,
15], Infomap [14], Collective Symmetric Nonnegative Matrix Factorization
(CSNMF) [30], Collective Projective Nonnegative Matrix Factorization (CPNMF)
[30], and Collective Symmetric Nonnegative Matrix Tri-factorization (CSNMTF)
[30].

####  VII-A3 Experiment 1

In this experiment, we generated two different types of networks, one where
the common communities are present across all layers and another where the
common communities are present in different subsets of layers. Fig. 3 shows
the results for the networks with 2 common communities across all layers for 3
(3(a)), 4 (3(b)), and 5 layers (3(c)), and for the networks with 3 common
communities across different subsets of layers for 3 (3(d)), 4 (3(e)), and 5
layers (3(f)). The results indicate that our method performs well for both
networks with common communities across all layers as well as for networks
with common communities that do not span all layers. Our method discovers the
complete structure of the network rather than forcing it to have a consensus
partition. Moreover, our method is robust to noise for larger values of $\mu$
compared to the other methods. We can also conclude that our algorithm
performs better when the common communities are across all layers (see Fig.
3(a), 3(b), 3(c)) than when the multiplex community structure is more complex
with common communities across subsets of layers (see Fig. 3(d), 3(e), 3(f) ),
but it still outperforms the rest of the methods. From Fig. 3 we can see that
GenLouvain performs well when $\mu$ is small, but its performance deteriorates
for $\mu$ values above 0.6. Another observation is that when the number of
layers is small, the NMF-based methods perform closer to GenLouvain but when
the number of layers increases, NMF algorithms perform worse. This is because
these NMF methods perform aggregation on either the adjacency or the community
indicator matrices. When there is more variation across layers, these methods
fail to capture this heterogeneity.

####  VII-A4 Experiment 2

In the second experiment, we evaluated the robustness of the algorithm against
variations in the common community structure by fixing $\mu=0.1$ and varying
the inter-layer dependency probability, $p_{1}$, i.e., the common communities
are allowed to vary across layers. The performance of all methods for a
5-layer network are reported in Fig. 4 based on the average NMI over 100
realizations of the network. As we can see in Fig. 4, our method still
outperforms the other seven methods when there is some variation in the common
community structure. This demonstrates that our method is robust to variations
of the common community structure across layers. However, our algorithm is
more sensitive to the drop in $p_{1}$ than the rest of the methods. This is
because when the common communities have a high variation our algorithm may
try to assign some of those nodes to the private communities.

Figure 4: 5-layer network generated with 6 different values of the interlayer
dependency probability $p_{1}$, with $\mu=0.1$, and $n=256$

####  VII-A5 Experiment 3

Another parameter in our model is the number of common communities $k_{c}$. In
this experiment, we fixed $n=256$, $\mu=0.3$, and the number of communities in
each layer, and varied $k_{c}$ from 1 to 7. When $k_{c}=7$, all communities
are common across layers and there are no private communities. As we can see
in Fig. 5, as $k_{c}$ is increased, the performance of all the other methods
improves, as expected, because these methods are designed to detect the common
community structure. When the communities are common across layers, most of
the methods, except Infomap, converge to the same NMI value. The performance
of MX-ONMTF is not affected by increasing $k_{c}$, and when all communities
are common it performs similarly to the other methods.

Figure 5: 5-layer network generated with different values of common
communities $k_{c}$ across layers.

####  VII-A6 Scalabilty Analysis

In this experiment, we evaluate the effect of network size on the run time of
the proposed algorithm. For this purpose, we fixed $\mu=0.3$, $L=5$, $k_{c}=3$
and varied $n$ from 32 to 8192. From Fig. 6, it can be seen that our method’s
run time is almost log-linear. This is comparable with all the other NMF based
methods. However, our as shown in the previous experiments, our method
performs better. Most of this time complexity is due to the multiplicative
update rule used in NMF-based algorithms and can be reduced using alternative
approaches as discussed in [50].

Figure 6: 5-layer network with 9 different values of $n$, $k_{c}=3$, and
$\mu=0.3$.

(a) Advice

(b) Friendship

(c) Co-work

Figure 7: Communities detected in Lazega Law Firm network across the three layers, advice, friendship, and co-work relationships. Red nodes are in the common community across the three layers. Method | Status | Gender | Office | Seniority | Age | Practice | Law School
---|---|---|---|---|---|---|---
GenLouvain | 0.0345 | 0.0307 | 0.5294 | 0.0807 | 0.0431 | 0.5468 | 0.0040
Aggregated Average | 0.0383 | 0.0197 | 0.5379 | 0.1307 | 0.0798 | 0.4411 | 0.0201
SC_ML | 0.0138 | 0.0259 | 0.0731 | 0.1225 | 0.0464 | 0.0249 | 0.0140
Infomap | 0.0179 | 0.0043 | 0.1668 | 0.2880 | 0.0083 | 0.0003 | 0.0093
CSNMF | 0.0418 | 0.0291 | 0.5732 | 0.1155 | 0.0736 | 0.4227 | 0.0172
CPNMF | 0.0081 | 0.0524 | 0.1139 | 0.0514 | 0.0291 | 0.0187 | 0.0221
CSNMTF | 0.0395 | 0.0217 | 0.0798 | 0.0795 | 0.0487 | 0.1335 | 0.0279
MX-ONMTF | 0.4752 | 0.4906 | 0.7386 | 0.4135 | 0.4203 | 0.6162 | 0.4226
TABLE III: NMI of the obtained community partition for each method and the
metadata available for the Lazega Law Firm Multiplex Network.

###  VII-B Real World Multiplex Networks

####  VII-B1 Lazega Law Firm Multiplex Social Network

Lazega Law Firm [51] is a multiplex social network with 71 nodes and three
layers representing Co-work, Friendship and Advice relationships between
partners and associates of a corporate law firm. This data set also includes
information about some attributes of each node such as status, gender, office
location, years with the firm, age, type of practice, and law school.

Applying MX-ONMTF to this network, we obtain one common community across all
layers composed of the nodes colored in red as well as private communities for
each layer, as shown in Fig. 7. This network does not have ground truth
community structure, but we can compute the NMI between the detected community
structure and each type of node attributes, i.e., metadata, to gain better
insight into the results and to be able to provide quantitative results [18].
For each of the attributes, the nodes are divided into communities based on
that particular attribute. For example, for the status, the network is divided
into two communities, partners and associates. For Age and Seniority, the
nodes were grouped into five-year bins. The community structure for each
attribute is used as ground truth to compute the NMI between each attribute
and the community structure detected by our method. The NMI values given in
Table III for the partition obtained by our method, suggest that office
location and type of practice (litigation or corporate) are highly correlated
with community membership across co-work, friendship and advice relationships.
We can also see that the partition detected by MX-ONMTF has greater NMI values
for each of the attributes. Therefore, our method detects a community
structure that takes all of the attributes into account instead of
partitioning with respect to just one attribute as the Aggregated Average
does.

####  VII-B2 C. Elegans Network

C. Elegans Network [52, 53] is a multiplex network with 279 nodes and 3 layers
representing different synaptic junctions (electric, chemical monadic, and
polyadic) of 279 neurons of the Caenorhabditis Elegans connectome. Information
about different attributes of the neurons in this dataset such as the group of
neuron they belong to (bodywall, mechanosensory, ring interneurons, head motor
neurons, etc.), the type of neuron (motor neurons, sensory neurons,
interneurons), and the color (blue, red, yellow, orange, etc.) is available.

Table IV shows the NMI values between the community structures detected by
each method and each of the three attributes available for this dataset. The
partition detected by MX-ONMTF has greater NMI values for each of the
attributes compared to the other four methods.

Method | Neuron Group | Neuron Type | Color
---|---|---|---
GenLouvain | 0.3756 | 0.1297 | 0.2362
Aggregated Average | 0.3839 | 0.1590 | 0.2977
SC-ML | 0.0185 | 0.0103 | 0.2690
Infomap | 0.2265 | 0.2345 | 0.2355
CSNMF | 0.1635 | 0.075 | 0.1211
CPNMF | 0.0854 | 0.0277 | 0.0628
CSNMTF | 0.1113 | 0.0402 | 0.0914
MX-ONMTF | 0.4074 | 0.4001 | 0.4593
TABLE IV: NMI of the obtained community partition for each method and the
metadata available for C. Elegans Network.

####  VII-B3 YeastLandscape Multiplex Network

Yeast Landscape is a multiplex genetic interaction network of a specie of
yeast, Saccharomyces Cerevisiae [54, 52]. This network has 4458 nodes and 4
layers representing the positive and negative interaction networks of genes in
Saccharomyces cerevisiae and positive and negative correlation based networks
in which genes with similar interaction profiles are connected to each other.
For this paper, we use the bioprocess annotations of the genes available on
the supplementary data file S6 of [54] as ground truth. We divided the genes
into 18 groups according to their primary bioprocess. There were 1580 genes in
this network without attributes.

Table V shows the NMI values between the community structures detected by each
method and the bioprocess of the genes. MX-ONMTF gives the highest NMI value
followed by the other NMF-based community detection methods.

Method | Bioprocess
---|---
GenLouvain | 0.0794
Aggregated Average | 0.1108
SC-ML | 0.1564
Infomap | 0.2987
CSNMF | 0.3553
CPNMF | 0.3559
CSNMTF | 0.3549
MX-ONMTF | 0.4123
TABLE V: NMI of the obtained community partition for each method with respect
to the metadata available for YeastLandscape Network.

###  VII-C Multiview Networks

In order to evaluate the performance of our method on networks where the
communities are common across all layers, we use two multiview data sets, UCI
Handwritten Digits111https://archive.ics.uci.edu/ml/datasets/Multiple+Features
[55] and Caltech [56].

The UCI Handwritten Digits data set consists of features of handwritten digits
from (0- 9) extracted from a collection of Dutch utility maps. There is a
total of 2000 patterns that have been digitized in binary images, 200 patterns
per digit. These digits are represented by six different feature sets: Fourier
coefficients of the character shapes, profile correlations, Karhunen-Loève
coefficients, pixel averages in 2 $\times$ 3 windows, Zernike moments, and
morphological features. Each layer of the multiplex network represents one of
the 6 features. The graphs are constructed using $k$-nearest neighbors graphs
with the nearest 50 neighbors and Euclidean distance.

Caltech-101 is a well-known object recognition dataset that consists of
pictures of objects belonging to 102 categories. There are about 40 to 800
images per category for a total of 9144 images. This dataset consists of 6
types of features extracted from each image. A multiplex network with 6 layers
representing each of the features, 102 classes, and 9144 nodes is constructed
from this dataset using k-nearest neighbors graphs with the nearest 50
neighbors. A smaller version of this dataset is also used in these
experiments, where only 20 objects are selected, resulting in a multiplex
network with 6 layers, 20 communities, and 2386 nodes.

In this case, as we have the true class assignment, we compute the NMI with
respect to this ground truth. As it can be seen in Table VI, our method
performs better than the rest of the methods for the three networks. This
indicates that even in cases where there are no private communities, our
method is successful at obtaining the consensus community structure, thus can
be used as an alternative to multiview clustering.

Method | Handwritten | Caltech-20 | Caltech-101
---|---|---|---
GenLouvain | 0.8791 | 0.5921 | 0.3406
Aggregated Average | 0.7957 | 0.5358 | 0.3941
SC-ML | 0.8435 | 0.6476 | 0.5016
Infomap | 0.5367 | 0.3876 | 0.2583
CSNMF | 0.4499 | 0.4250 | 0.3842
CPNMF | 0.4421 | 0.4208 | 0.3816
CSNMTF | 0.4478 | 0.4264 | 0.3831
MX-ONMTF | 0.9432 | 0.6861 | 0.5660
TABLE VI: NMI of the obtained community partition for multiview networks.

##  VIII Conclusions

In this paper, we proposed a multiplex community detection method based on
ONMTF. The proposed method, MX-ONMTF, is able to detect both common and
private communities across layers, allowing us to differentiate between the
topologies across layers. The proposed algorithm is based on multiplicative
update rules and a proof of convergence is provided. A new approach based on
the eigengap criterion is introduced for determining the number of
communities. Results for both synthetic and real-world networks show that our
method performs better than existing community detection methods for multiplex
networks as it is able to handle the heterogeneity of the network topology
across layers. Moreover, experiments on multiview networks show that our
method also performs well in cases where a consensus community structure is
needed.

##  Appendix A Auxiliary Function Proof

######  Proposition 1.

The following function, $Z(h,h_{ij}^{t})$,

$\displaystyle\begin{split}&Z(h,h_{ij}^{t})=\mathcal{L}(h^{t}_{ij})+3\mathcal{L^{\prime}}(h^{t}_{ij})(h-h^{t}_{ij})+\\\
&\frac{3}{2}\frac{\sum_{l=1}^{L}(4\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l}+4\mathbf{H}^{t}{\mathbf{H}^{t}}^{\top}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}(h-h^{t}_{ij})^{2}\end{split}$

is an auxiliary function of $\mathcal{L}(H)$,

$\displaystyle\begin{split}\mathcal{L}(h)=\mathcal{L}(h^{t}_{ij})+\mathcal{L^{\prime}}(h^{t}_{ij})(h-h^{t}_{ij})+\frac{1}{2}\mathcal{L^{\prime\prime}}(h_{ij}^{t})(h-h^{t}_{ij})^{2}.\end{split}$

###### Proof.

First, when $h=h_{ij}^{t}$ the equality $Z(h,h)=\mathcal{L}(h)$ holds. Now, we
need to show that $Z(h,h_{ij}^{t})\geq\mathcal{L}(h)$.

It can be seen that the first and second terms of $Z(h,h_{ij}^{t})$ are
greater than the first and second terms in $\mathcal{L}(h)$. Therefore, it
suffices to show that
$\frac{3\sum_{l=1}^{L}(4\mathbf{H}_{l}\mathbf{G}_{l}{\mathbf{H}_{l}}^{\top}\mathbf{H}^{t}\mathbf{S}_{l}+4\mathbf{H}^{t}{\mathbf{H}^{t}}^{\top}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}\geq\mathcal{L^{\prime\prime}}(h_{ij}^{t})$.

It can be shown that

$\begin{split}\frac{(\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}=\frac{\sum_{p,q}(\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top})_{ip}h_{pq}^{qj}}{h_{ij}^{t}}\\\
\geq(\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top})_{ii}\mathbf{S}_{jj},\end{split}$

$\begin{split}\frac{(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}=\frac{\sum_{p}h_{ip}^{t}({\mathbf{H}^{t}}^{\top}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{pj}}{h_{ij}^{t}}\\\
\geq{({\mathbf{H}^{t}}^{\top}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{jj}},\end{split}$

$\displaystyle\frac{(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}=\frac{\sum_{p,q}h_{ip}^{t}h_{qp}^{t}(\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{qj}}{h_{ij}^{t}}$
$\displaystyle\geq
h_{ij}^{t}(\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij},$

$\displaystyle\frac{(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}=\frac{\sum_{m,b}(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l})_{im}h_{mb}^{bj}}{h_{ij}^{t}}$
$\displaystyle\geq(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l})_{ii}\mathbf{S}_{l_{jj}}.$

Therefore,

$\displaystyle\begin{split}3\frac{(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}&\geq({\mathbf{H}^{t}}^{\top}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}+h_{ij}^{t}(\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}\\\
+&(\mathbf{H}^{t}{\mathbf{H}}^{t^{\top}}\mathbf{A}_{l})_{ii}\mathbf{S}_{l_{jj}},\end{split}$

and thus
$\frac{3\sum_{l=1}^{L}(4\mathbf{H}_{l}\mathbf{G}_{l}\mathbf{H}_{l}^{\top}\mathbf{H}^{t}\mathbf{S}_{l}+4\mathbf{H}^{t}\mathbf{H}^{t^{\top}}\mathbf{A}_{l}\mathbf{H}^{t}\mathbf{S}_{l})_{ij}}{h_{ij}^{t}}\geq\mathcal{L^{\prime\prime}}(h_{ij}^{t})$.
Therefore, Eq. (11) is an auxiliary function of $\mathcal{L}(h)$. ∎

## Acknowledgment

This work was supported in part by National Science Foundation under Grants
CCF-2006800.

## References

  * [1] A.-L. Barabási, “Network science,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 371, no. 1987, p. 20120375, 2013.
  * [2] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, “Multilayer networks,” _Journal of complex networks_ , vol. 2, no. 3, pp. 203–271, 2014.
  * [3] S. E. Smith-Aguilar, F. Aureli, L. Busia, C. Schaffner, and G. Ramos-Fernández, “Using multiplex networks to capture the multidimensional nature of social structure,” _Primates_ , vol. 60, no. 3, pp. 277–295, 2019.
  * [4] A. Aleta, S. Meloni, and Y. Moreno, “A multilayer perspective for the analysis of urban transportation systems,” _Scientific reports_ , vol. 7, no. 1, pp. 1–9, 2017.
  * [5] S. Fortunato and D. Hric, “Community detection in networks: A user guide,” _Physics reports_ , vol. 659, pp. 1–44, 2016.
  * [6] M. Magnani, O. Hanteer, R. Interdonato, L. Rossi, and A. Tagarelli, “Community detection in multiplex networks,” _ACM Computing Surveys (CSUR)_ , vol. 54, no. 3, pp. 1–35, 2021.
  * [7] M. Berlingerio, M. Coscia, and F. Giannotti, “Finding and characterizing communities in multidimensional networks,” in _2011 International Conference on advances in social networks analysis and mining_. IEEE, 2011, pp. 490–494.
  * [8] P.-Y. Chen and A. O. Hero, “Multilayer spectral graph clustering via convex layer aggregation: Theory and algorithms,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 3, no. 3, pp. 553–567, 2017.
  * [9] D. Taylor, R. S. Caceres, and P. J. Mucha, “Super-resolution community detection for layer-aggregated multilayer networks,” _Physical Review X_ , vol. 7, no. 3, p. 031056, 2017.
  * [10] M. Berlingerio, F. Pinelli, and F. Calabrese, “Abacus: frequent pattern mining-based community discovery in multidimensional networks,” _Data Mining and Knowledge Discovery_ , vol. 27, no. 3, pp. 294–320, 2013.
  * [11] L. Tang, X. Wang, and H. Liu, “Community detection via heterogeneous interaction analysis,” _Data mining and knowledge discovery_ , vol. 25, no. 1, pp. 1–33, 2012.
  * [12] X. Dong, P. Frossard, P. Vandergheynst, and N. Nefedov, “Clustering on multi-layer graphs via subspace analysis on grassmann manifolds,” _IEEE Transactions on signal processing_ , vol. 62, no. 4, pp. 905–918, 2013.
  * [13] Z. Kuncheva and G. Montana, “Community detection in multiplex networks using locally adaptive random walks,” in _Proceedings of the 2015 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining 2015_ , 2015, pp. 1308–1315.
  * [14] M. De Domenico, A. Lancichinetti, A. Arenas, and M. Rosvall, “Identifying modular flows on multilayer networks reveals highly overlapping organization in interconnected systems,” _Physical Review X_ , vol. 5, no. 1, p. 011027, 2015.
  * [15] P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, “Community structure in time-dependent, multiscale, and multiplex networks,” _science_ , vol. 328, no. 5980, pp. 876–878, 2010.
  * [16] A. Amelio and C. Pizzuti, “Community detection in multidimensional networks,” in _2014 IEEE 26th International Conference on Tools with Artificial Intelligence_. IEEE, 2014, pp. 352–359.
  * [17] G. Zhu and K. Li, “A unified model for community detection of multiplex networks,” in _International Conference on Web Information Systems Engineering_. Springer, 2014, pp. 31–46.
  * [18] A. R. Pamfil, S. D. Howison, R. Lambiotte, and M. A. Porter, “Relating modularity maximization and stochastic block models in multilayer networks,” _SIAM Journal on Mathematics of Data Science_ , vol. 1, no. 4, pp. 667–698, 2019.
  * [19] M. Ortiz-Bouza and S. Aviyente, “Orthogonal nonnegative matrix tri-factorization for community detection in multiplex networks,” in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 5987–5991.
  * [20] H. T. Ali, S. Liu, Y. Yilmaz, R. Couillet, I. Rajapakse, and A. Hero, “Latent heterogeneous multilayer community detection,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 8142–8146.
  * [21] C. De Bacco, E. A. Power, D. B. Larremore, and C. Moore, “Community detection, link prediction, and layer interdependence in multilayer networks,” _Physical Review E_ , vol. 95, no. 4, p. 042317, 2017.
  * [22] O. Boutemine and M. Bouguessa, “Mining community structures in multidimensional networks,” _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , vol. 11, no. 4, pp. 1–36, 2017.
  * [23] F. Karimi, S. Lotfi, and H. Izadkhah, “Multiplex community detection in complex networks using an evolutionary approach,” _Expert Systems with Applications_ , vol. 146, p. 113184, 2020.
  * [24] C. Chen, M. K. Ng, and S. Zhang, “Block spectral clustering methods for multiple graphs,” _Numerical Linear Algebra with Applications_ , vol. 24, no. 1, p. e2075, 2017.
  * [25] C. Ding, X. He, and H. D. Simon, “On the equivalence of nonnegative matrix factorization and spectral clustering,” in _Proceedings of the 2005 SIAM international conference on data mining_ , 2005, pp. 606–610.
  * [26] H. Wang, F. Nie, H. Huang, and C. Ding, “Nonnegative matrix tri-factorization based high-order co-clustering and its fast implementation,” in _2011 IEEE 11th international conference on data mining_. IEEE, 2011, pp. 774–783.
  * [27] S. Mankad and G. Michailidis, “Structural and functional discovery in dynamic networks with non-negative matrix factorization,” _Physical Review E_ , vol. 88, no. 4, p. 042812, 2013.
  * [28] B.-J. Sun, H. Shen, J. Gao, W. Ouyang, and X. Cheng, “A non-negative symmetric encoder-decoder approach for community detection,” in _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management_ , 2017, pp. 597–606.
  * [29] X. Ma, D. Dong, and Q. Wang, “Community detection in multi-layer networks using joint nonnegative matrix factorization,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 31, no. 2, pp. 273–286, 2018.
  * [30] V. Gligorijević, Y. Panagakis, and S. Zafeiriou, “Non-negative matrix factorizations for multiplex network analysis,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 41, no. 4, pp. 928–940, 2018\.
  * [31] H. T. Nguyen, T. N. Dinh, and T. Vu, “Community detection in multiplex social networks,” in _2015 IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS)_. IEEE, 2015, pp. 654–659.
  * [32] U. Von Luxburg, “A tutorial on spectral clustering,” _Statistics and computing_ , vol. 17, no. 4, pp. 395–416, 2007.
  * [33] E. Cozzo, M. Kivelä, M. De Domenico, A. Solé-Ribalta, A. Arenas, S. Gómez, M. A. Porter, and Y. Moreno, “Structure of triadic relations in multiplex networks,” _New Journal of Physics_ , vol. 17, no. 7, p. 073029, 2015.
  * [34] W. Wu, S. Kwong, Y. Zhou, Y. Jia, and W. Gao, “Nonnegative matrix factorization with mixed hypergraph regularization for community detection,” _Information Sciences_ , vol. 435, pp. 263–281, 2018.
  * [35] H. Lu, X. Sang, Q. Zhao, and J. Lu, “Community detection algorithm based on nonnegative matrix factorization and pairwise constraints,” _Physica A: Statistical Mechanics and its Applications_ , vol. 545, p. 123491, 2020.
  * [36] C. Ding, T. Li, W. Peng, and H. Park, “Orthogonal nonnegative matrix t-factorizations for clustering,” in _Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2006, pp. 126–135.
  * [37] D. Lee and H. Seung, “Algorithms for non-negative matrix factorization,” vol. 13, pp. 535–541, 2001.
  * [38] J. Yoo and S. Choi, “Orthogonal nonnegative matrix tri-factorization for co-clustering: Multiplicative updates on stiefel manifolds,” _Information processing & management_, vol. 46, no. 5, pp. 559–570, 2010\.
  * [39] T. Li and C.-c. Ding, “Nonnegative matrix factorizations for clustering: A survey,” in _Data Clustering_. Chapman and Hall/CRC, 2018, pp. 149–176.
  * [40] X. Luo, Z. Liu, L. Jin, Y. Zhou, and M. Zhou, “Symmetric nonnegative matrix factorization-based community detection models and their convergence analysis,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 33, no. 3, pp. 1203–1215, 2021.
  * [41] L. Danon, A. Diaz-Guilera, J. Duch, and A. Arenas, “Comparing community structure identification,” _Journal of statistical mechanics: Theory and experiment_ , vol. 2005, no. 09, p. P09008, 2005.
  * [42] Z. Li, S. Zhang, R.-S. Wang, X.-S. Zhang, and L. Chen, “Quantitative function for community detection,” _Physical review E_ , vol. 77, no. 3, p. 036109, 2008.
  * [43] S. Pramanik, R. Tackx, A. Navelkar, J.-L. Guillaume, and B. Mitra, “Discovering community structure in multilayer networks,” in _2017 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_. IEEE, 2017, pp. 611–620.
  * [44] F. Liu, D. Choi, L. Xie, and K. Roeder, “Global spectral clustering in dynamic networks,” _Proceedings of the National Academy of Sciences_ , vol. 115, no. 5, pp. 927–932, 2018.
  * [45] S. Paul and Y. Chen, “Spectral and matrix factorization methods for consistent community detection in multi-layer networks,” _The Annals of Statistics_ , vol. 48, no. 1, pp. 230–250, 2020.
  * [46] M. Bazzi, L. Jeub, A. Arenas, S. D. Howison, and M. A. Porter, “Generative benchmark models for mesoscale structure in multilayer networks,” _arXiv preprint arXiv:1608.06196_ , p. 20, 2016.
  * [47] L. Jeub and M. Bazzi, “A generative model for mesoscale structure in multilayer networks implemented in matlab,” 2016. [Online]. Available: https://github.com/MultilayerGM/MultilayerGM-MATLAB
  * [48] B. Karrer and M. E. Newman, “Stochastic blockmodels and community structure in networks,” _Physical review E_ , vol. 83, no. 1, p. 016107, 2011.
  * [49] I. S. Jutla, L. G. Jeub, P. J. Mucha _et al._ , “A generalized louvain method for community detection implemented in matlab.” [Online]. Available: https://github.com/GenLouvain/GenLouvain (2011-2019)
  * [50] A. Čopar, B. Zupan, and M. Zitnik, “Fast optimization of non-negative matrix tri-factorization,” _PloS one_ , vol. 14, no. 6, p. e0217994, 2019\.
  * [51] E. Lazega _et al._ , _The collegial phenomenon: The social mechanisms of cooperation among peers in a corporate law partnership_. Oxford University Press on Demand, 2001.
  * [52] M. De Domenico, M. A. Porter, and A. Arenas, “Muxviz: a tool for multilayer analysis and visualization of networks,” _Journal of Complex Networks_ , vol. 3, no. 2, pp. 159–176, 2015.
  * [53] B. L. Chen, D. H. Hall, and D. B. Chklovskii, “Wiring optimization can relate neuronal structure and function,” _Proceedings of the National Academy of Sciences_ , vol. 103, no. 12, pp. 4723–4728, 2006.
  * [54] M. Costanzo, A. Baryshnikova, J. Bellay, Y. Kim, E. D. Spear, C. S. Sevier, H. Ding, J. L. Koh, K. Toufighi, S. Mostafavi _et al._ , “The genetic landscape of a cell,” _science_ , vol. 327, no. 5964, pp. 425–431, 2010\.
  * [55] D. Dua and C. Graff, “Uci machine learning repository,” 2017. [Online]. Available: http://archive.ics.uci.edu/ml
  * [56] Li F. F., Andreeto M., Ranzato M., and Perona P., “Caltech 101 (1.0) [Data set]. CaltechDATA.” 2022. [Online]. Available: https://doi.org/10.22002/D1.20086


Rumen Bogdanovski, Nobody Else 11institutetext:  Institute of Astronomy and
NAO, Bulgarian Academy of Sciences, BG-1784, Sofia 22institutetext: Shumen
University, 115, Universitetska St, BG-9700, Shumen
22email<EMAIL_ADDRESS>Submitted on xx.xx.xxxx; Accepted on xx.xx.xxxx

# New Results from Long-time Photometric Study
of UX Orionis Star GM Cephei

Asen Mutafov 11 Evgeni Semkov 11 Stoyanka Peneva 11 Sunay Ibryamov 22

###### Abstract

New results from long-term optical photometric observations of the pre-main
sequence star GM Cep from UX Orionis type are reported. During ongoing
photometric monitoring of the GM Cep four deep minimums in brightness are
observed. The collected multicolour photometric data shows the typical of UXor
variables colour reversal during the minimums in brightness. Recent $BVRI$
photometric observations of GM Cep have been collected from November 2014 to
October 2020.

###### keywords:

stars, stellar evolution, UX Orionis stars

## 1\. Introduction

The pre-main sequence star (PMS) GM Cep is situated in the field of the young
open cluster Trumpler 37 ($\sim$4 Myr old and at a distance of 870 pc,
Contreras et al. [3]). It is likely to be a member of the cluster, Marschall
and van Altena [5], Sicilia-Aguilar et al. [8].

According to the study of Sicilia-Aguilar et al. [9] GM Cep has a solar-like
mass ($\it M$ $\sim$ $2.1M_{\odot}$), with radius that ranges between 3 and 6
$R_{\odot}$. It is from G7V-K0V spectral type and with a strong IR excesses
that can be explained by the presence of a very luminous and massive
circumstellar disk. In the study is found out that the H$\alpha$ emission line
in the spectrum of GM Cep has a strong P Cyg profile. It is also observed that
the equivalent width of the line varies significantly from 6$\AA$ to 19$\AA$
and that the variable accretion rate is up to $\sim$ 10-6 $M_{\odot}$/year.

A long-term photometric study of GM Cep was made by Xiao et al. [10] in the
period of several decades. The long-term B and V light curves of the star are
made by using the photographic plate archives from Harvard College Observatory
and from Sonneberg Observatory. The results suggest that the light curves of
GM Cep seem to be dominated by dips superposed on the quiescence state and
that there is a lack of fast rises in brightness typical of EXor variables. In
their study, Xiao et al. [10] have not found evidence for periodicity of
observed dips in brightness.

The results from the $BVRI$ photometric observations of GM Cep collected in
the period 2008 June - 2014 August and reported in our previous studies
(Semkov and Peneva[6] and Semkov et al. [7]) show very strong photometric
variability of the star. In this period we have registered five deep minimums
in brightness in the light curve of GM Cep. On the basis of these observations
back then we concluded that the variability of GM Cep is dominated by fading
events rather than by bursting events. The collected multicolour photometric
data shows the effect of a colour reversal at the deep minimum of brightness,
which is evidence of variable extinction from the circumstellar environment,
typical of UXor variables.

In the study of Chen et al. [1] while carrying out an intensive $BVR$
photometric monitoring of GM Cep during the period 2009-2011 they confirm the
UXor nature of its variability and suggest an early stage of planetesimal
formation in the star environment. A periodicity of about 300 days of the
observed deep declines in brightness is suggested by Chen and Hu [2].

The multicolour observations give us the opportunity to clarify the mechanism
of the brightness variations.

## 2\. Observations

The CCD observations of GM Cep cover the period from June 2008 to October
2020. They were performed in two observatories with four telescopes: Rozhen
National Astronomical Observatory (Bulgaria) with its 2-metre RCC, 50/70-cm
Schmidt and 60-cm Cassegrain telescopes along with Skinakas Observatory of the
University of Crete (Greece) with the 1.3-m Ritchey-Crétien telescope. Five
different types of CCD cameras were used during the observations. Their
technical characteristics and optical specifications are given in Table 1.

We used the published in the work of Semkov and Peneva [6] fifteen stars in
the field around GM Cep for reference. A standard set of Johnson-Cousins’
filters were used for all the frames that were taken. Twilight flat fields in
each filter were obtained each clear evening. All frames obtained with the
ANDOR and Vers Array cameras are bias subtracted and flat fielded. CCD frames
obtained with the FLI PL16803 and FLI PL09000 cameras are dark subtracted and
flat fielded. Using IDL DAOPHOT routines an aperture photometry was performed.
In order to obtain comparable results with our previous obserations we used
the same aperture to analyse all the data. It was chosen as 6 arcsec in
radius, while the background annulus was from 10 to 15 arcsec.

Table 1: CCD cameras and optical specifications Telescope | CCD Camera Type | Size | Field | Pixel | Scale | RON | Gain
---|---|---|---|---|---|---|---
|  |  |  | size |  |  |
|  | (px) | (arcmin) |  ($\mu$m) | (”/px) |  ($e^{-}$rms) |  ($e^{-}$/ADU)
2m RCC | Vers Array 1300B | 1340x1300 | 5.8x5.6 | 20.0 | 0.26 | 2.00 | 1.0
2m RCC | ANDOR iKon-L | 2048x2048 | 6.0x6.0 | 13.5 | 0.17 | 6.90 | 1.1
Schmidt | FLI PL16803 | 4096x4096 | 73.8x73.8 | 9.0 | 1.08 | 9.00 | 1.0
60cm Cass | FLI PL9000 | 3056x3056 | 16.8x16.8 | 12.0 | 0.33 | 8.50 | 1.0
1.3m RC | ANDOR DZ436-BV | 2048x2048 | 9.6x9.6 | 13.5 | 0.28 | 8.14 | 2.7

## 3\. Results

In the current paper the data from the multicolour photometric observations of
GM Cep is presented for the period from August 2014 to October 2020. They are
a continuation of the observations that were begun by us in June 2008.
Previous data is published in the work of Semkov and Peneva [6] for the period
June 2008 - February 2011 and in the work of Semkov et al. [7] for the period
April 2011 - August 2014.

The new data in the period from November 2014 to October 2020 is shown in
Table 2. The columns provide the Julian date (JD) of observation, $BVRI$
magnitudes, and the telescope used. In the column Telescope the abbreviation
2-m denotes the 2-m Ritchey-Chrétien-Coudé, Schmidt - the 50/70-cm Schmidt and
the 1.3-m Ritchey-Crétien telescope. The values of the instrumental errors are
in the range 0.m01-0.m05 (for B), 0.m01-0.m03 (for V) and 0.m01-0.m02 (for R
and I) (Semkov and Peneva [6]).

As presented by the graphics in Figure 1 the new photometric data shows
continued strong brightness variability of GM Cep. The same brightness
variability was also registered in the previous studies Sicilia-Aguilar et al.
[9], Xiao et al. [10], Semkov and Peneva [6], Chen et al. [1], Semkov et al.
[7], Huang et al. [4].

In the time scale of days and months outside the deep minimums GM Cep also
shows significant brightness variations. The summarized results of over 12
years period of observations show very strong photometric variability. We have
registered new four deep minimums in brightness in the light curve of GM Cep:
August 2015, January 2017, November 2017, August 2020.

In Figure 2 are shown respectively colour-magnitude diagrams $V$/$B-V$,
$V$/$V-R$ and $V$/$V-I$. The collected multicolour photometric data shows the
typical of UXor colour reversal during the minimums in brightness. This
coincides with the model of blurring of the dust-like material. The observed
reverse of the colour is caused by the diffused light from small dust grains.
This is a typical characteristic of PMS stars from the Uxor type. From visual
inspection for every of the colour diagrams such point of reversal is observed
at different star brightness: in the $V$/$B-V$ diagram, the point of reversal
is observed at $V$ about 14.0 mag, in the $V$/$V-R$ diagram at $V$ about 14.5
mag and in the $V$/$V-I$ diagram at $V$ about 14.6 mag.

Usually, when in the line of sight there are clusters of dust the star becomes
redder. But during a maximum eclipse the blue part of the diffused light in
the observed light becomes significant and the star colour becomes bluer.

Figure 1: BVRI light curves of GM Cep for the whole period of our photometric monitoring (2008 - 2020). Figure 2: The colour-magnitude diagrams $V$/$B-V$, $V$/$V-R$ and $V$/$V-I$ of GM Cep in the period of observations June 2008 - October 2020. Table 2: $BVRI$ photometric observations of GM Cep JD (24…) | $\it B$ | $\it V$ | $\it R$ | $\it I$ | Tel | JD (24…) | $\it B$ | $\it V$ | $\it R$ | $\it I$ | Tel
---|---|---|---|---|---|---|---|---|---|---|---
56988.217 | 14.51 | 13.26 | 12.48 | 11.57 | Sch | 57969.371 | 14.46 | 13.17 | 12.39 | 11.55 | Sch
57005.259 | - | 13.15 | 12.33 | 11.45 | Sch | 58011.342 | 14.81 | 13.51 | 12.65 | 11.71 | Sch
57006.315 | 14.34 | 13.14 | 12.32 | 11.45 | Sch | 58012.352 | 14.93 | 13.58 | 12.72 | 11.74 | Sch
57016.251 | 14.41 | 13.16 | 12.29 | 11.49 | 2-m | 58013.348 | 14.88 | 13.55 | 12.69 | 11.75 | Sch
57017.202 | 14.40 | 13.10 | 12.28 | 11.51 | 2-m | 58039.325 | - | - | 12.51 | 11.61 | Sch
57074.579 | 14.60 | 13.32 | 12.52 | 11.59 | Sch | 58043.304 | 14.59 | 13.27 | 12.48 | 11.62 | Sch
57136.511 | 14.71 | 13.41 | 12.60 | 11.66 | Sch | 58080.294 | 16.27 | 14.81 | 13.83 | 12.66 | Sch
57138.450 | 14.54 | 13.27 | 12.47 | 11.56 | Sch | 58081.311 | 16.39 | 14.95 | 13.98 | 12.85 | Sch
57162.406 | 14.47 | 13.23 | 12.42 | 11.52 | Sch | 58109.346 | 15.61 | 14.18 | 13.29 | 12.20 | Sch
57164.479 | 14.58 | 13.31 | 12.47 | 11.59 | Sch | 58113.284 | 15.50 | 14.08 | 13.16 | 12.09 | Sch
57167.433 | 14.52 | 13.21 | 12.41 | 11.57 | 2-m | 58114.287 | 15.63 | 14.18 | 13.27 | 12.18 | Sch
57186.480 | 14.66 | 13.40 | 12.55 | 11.64 | Sch | 58217.551 | 15.08 | 13.79 | 12.96 | 11.97 | Sch
57187.498 | 14.68 | 13.39 | 12.54 | 11.66 | 2-m | 58218.535 | 15.09 | 13.79 | 12.95 | 11.95 | Sch
57190.410 | 14.48 | 13.16 | 12.37 | 11.53 | 2-m | 58220.493 | 15.00 | 13.71 | 12.89 | 11.91 | Sch
57220.402 | 14.30 | 13.02 | 12.21 | 11.40 | Sch | 58278.436 | 15.19 | 13.84 | 12.94 | 11.91 | Sch
57221.459 | 14.22 | 13.05 | 12.25 | 11.44 | Sch | 58312.365 | 15.27 | 13.90 | 13.02 | 12.00 | Sch
57223.437 | 14.37 | 13.08 | 12.29 | 11.50 | 2-m | 58316.323 | 15.09 | 13.74 | 12.85 | 11.85 | Sch
57246.395 | 14.83 | 13.48 | 12.62 | 11.74 | 1.3-m | 58340.393 | 14.62 | 13.33 | 12.51 | 11.61 | Sch
57247.411 | 14.71 | 13.37 | 12.52 | 11.66 | 1.3-m | 58342.379 | 14.58 | 13.28 | 12.46 | 11.58 | Sch
57259.363 | 15.86 | 14.67 | 13.71 | 12.59 | Sch | 58343.391 | 14.58 | 13.29 | 12.47 | 11.59 | Sch
57260.404 | 15.93 | 14.65 | 13.81 | 12.81 | Sch | 58344.388 | 14.55 | 13.26 | 12.45 | 11.57 | Sch
57269.401 | 15.45 | 14.03 | 13.05 | 11.94 | Sch | 58346.344 | 14.62 | 13.36 | 12.53 | 11.64 | 2-m
57271.468 | 15.04 | 13.57 | 12.63 | 11.71 | 2-m | 58363.360 | 14.94 | 13.60 | 12.74 | 11.82 | Sch
57330.255 | 15.20 | 13.87 | 12.93 | 11.89 | Sch | 58364.323 | 14.83 | 13.50 | 12.66 | 11.73 | Sch
57331.274 | 15.26 | 13.93 | 13.00 | 11.97 | Sch | 58365.533 | 14.90 | 13.56 | 12.69 | 11.80 | 2-m
57332.264 | 15.34 | 14.00 | 13.08 | 12.05 | Sch | 58409.220 | 14.60 | 13.29 | 12.49 | 11.64 | Sch
57333.265 | 15.39 | 14.04 | 13.12 | 12.08 | Sch | 58428.174 | 15.34 | 13.90 | 12.97 | 11.95 | Sch
57334.244 | 15.41 | 14.08 | 13.18 | 12.11 | Sch | 58435.295 | 15.49 | 14.04 | 13.09 | 12.04 | Sch
57369.257 | 14.63 | 13.29 | 12.47 | 11.61 | 2-m | 58492.225 | 15.16 | 13.75 | 12.84 | 11.84 | Sch
57370.229 | 14.62 | 13.23 | 12.45 | 11.59 | 2-m | 58496.201 | 15.24 | 13.87 | 12.94 | 11.92 | Sch
57371.224 | 14.59 | 13.27 | 12.44 | 11.61 | 2-m | 58547.625 | 14.68 | 13.36 | 12.54 | 11.65 | Sch
57372.241 | 14.63 | 13.34 | 12.47 | 11.58 | Sch | 58603.481 | 15.14 | 13.74 | 12.83 | 11.84 | Sch
57374.299 | 14.73 | 13.42 | 12.51 | 11.58 | Sch | 58604.532 | 15.07 | 13.68 | 12.77 | 11.79 | Sch
57425.221 | 15.18 | 13.82 | 12.94 | 11.94 | Sch | 58665.455 | 14.25 | 12.98 | 12.22 | 11.41 | Sch
57426.226 | 15.22 | 13.88 | 12.97 | 11.99 | Sch | 58666.459 | 14.38 | 13.10 | 12.31 | 11.47 | Sch
57483.462 | 14.47 | 13.19 | 12.39 | 11.52 | 2-m | 58667.494 | 14.40 | 13.12 | 12.33 | 11.49 | Sch
57484.480 | 14.44 | 13.20 | 12.38 | 11.52 | 2-m | 58690.370 | 14.38 | 13.08 | 12.31 | 11.50 | 2-m
57522.460 | 14.69 | 13.39 | 12.58 | 11.63 | Sch | 58691.450 | 14.39 | 13.06 | 12.31 | 11.51 | 2-m
57523.448 | 14.70 | 13.39 | 12.57 | 11.63 | Sch | 58692.389 | 14.44 | 13.06 | 12.31 | 11.54 | 2-m
57540.447 | 14.89 | 13.43 | 12.55 | 11.65 | 2-m | 58704.371 | 14.72 | 13.40 | 12.56 | 11.66 | Sch
57565.485 | 14.92 | 13.56 | 12.69 | 11.67 | Sch | 58705.404 | 14.60 | 13.29 | 12.46 | 11.60 | Sch
57581.424 | 14.75 | 13.41 | 12.58 | 11.64 | Sch | 58706.402 | 14.53 | 13.23 | 12.41 | 11.55 | Sch
57582.459 | 14.67 | 13.38 | 12.56 | 11.62 | Sch | 58707.476 | 14.53 | 13.21 | 12.40 | 11.54 | Sch
57583.434 | 14.74 | 13.41 | 12.57 | 11.63 | Sch | 58726.384 | 14.44 | 13.14 | 12.44 | 11.55 | 2-m
57603.384 | 14.88 | 13.47 | 12.63 | 11.74 | 2-m | 58727.428 | 14.38 | 13.07 | 12.28 | 11.51 | 2-m
57605.412 | 14.88 | 13.54 | 12.72 | 11.74 | Sch | 58728.465 | 14.33 | 13.00 | 12.26 | 11.46 | 2-m
57607.392 | 14.72 | 13.40 | 12.60 | 11.65 | Sch | 58729.395 | 14.32 | 13.02 | 12.27 | 11.47 | 2-m
57664.304 | 14.57 | 13.28 | 12.52 | 11.63 | Sch | 58730.464 | 14.38 | 13.09 | 12.30 | 11.46 | Sch
57698.280 | 14.31 | 13.07 | 12.32 | 11.47 | Sch | 58758.414 | 14.25 | 12.97 | 12.20 | 11.38 | Sch
57714.306 | 14.48 | 13.16 | 12.38 | 11.57 | 2-m | 58759.461 | 14.29 | 13.01 | 12.23 | 11.41 | Sch
57715.272 | 14.39 | 13.12 | 12.35 | 11.53 | 2-m | 58864.243 | 14.22 | 12.97 | 12.20 | 11.39 | Sch
57716.286 | 14.38 | 13.10 | 12.32 | 11.50 | 2-m | 58865.243 | 14.20 | 12.96 | 12.18 | 11.39 | Sch
57755.231 | 15.68 | 14.29 | 13.39 | 12.29 | Sch | 58869.221 | 14.33 | 13.05 | 12.26 | 11.48 | 2-m
57756.243 | 15.64 | 14.26 | 13.34 | 12.24 | Sch | 58870.238 | 14.34 | 13.07 | 12.28 | 11.46 | Sch
57781.221 | 15.16 | 13.81 | 12.93 | 11.87 | Sch | 58993.457 | 14.49 | 13.18 | 12.35 | 11.49 | Sch
57782.230 | 15.19 | 13.82 | 12.85 | 11.90 | 2-m | 59040.398 | 15.34 | 14.00 | 13.07 | 12.07 | Sch
57784.227 | - | 13.75 | 12.81 | 11.86 | 2-m | 59041.413 | 15.36 | 14.02 | 13.09 | 12.09 | Sch
57785.251 | 15.17 | 13.75 | 12.83 | 11.89 | 2-m | 59042.403 | 15.33 | 14.02 | 13.09 | 12.10 | Sch
57786.222 | 15.20 | 13.77 | 12.86 | 11.91 | 2-m | 59059.412 | 15.37 | 14.10 | 13.21 | 12.22 | Sch
57800.205 | 15.00 | 13.67 | 12.81 | 11.80 | Sch | 59060.436 | 15.39 | 14.11 | 13.21 | 12.21 | Sch
57801.208 | 14.98 | 13.64 | 12.78 | 11.79 | Sch | 59075.371 | 15.61 | 14.25 | 13.31 | 12.28 | 2-m
57817.531 | 14.54 | 13.28 | 12.51 | 11.57 | Sch | 59085.314 | 15.39 | 14.03 | 13.10 | 12.10 | Sch
57845.516 | 14.35 | 13.15 | 12.39 | 11.51 | Sch | 59101.456 | 15.27 | 13.91 | 12.98 | 11.98 | Sch
57846.547 | 14.57 | 13.34 | 12.57 | 11.68 | Sch | 59102.383 | 15.27 | 13.87 | 12.94 | 11.95 | 2-m
57892.498 | 14.76 | 13.50 | 12.71 | 11.79 | Sch | 59103.464 | 15.32 | 13.91 | 12.99 | 12.00 | 2-m
57893.573 | 14.83 | 13.53 | 12.72 | 11.87 | 2-m | 59105.338 | 15.23 | 13.89 | 12.96 | 11.97 | Sch
57904.451 | 14.56 | 13.30 | 12.50 | 11.65 | Sch | 59108.400 | 14.78 | 13.55 | 12.68 | 11.76 | Sch
57967.438 | - | - | 12.50 | 11.60 | Sch | 59109.350 | 14.84 | 13.52 | 12.64 | 11.73 | Sch
57968.369 | 14.54 | 13.27 | 12.49 | 11.61 | Sch | 59136.280 | 14.76 | 13.44 | 12.58 | 11.67 | Sch

## 4\. Conclusion

In the time scale of days and months outside the deep minimums GM Cep shows
significant brightness variations.

The last data confirms again the presence of an “blueing effect” during the
minimum of brightness and is an independent evidence that the variability of
GM Cep is dominated by the variable extinction. The collected multicolour
photometric data shows the typical of UXor variables colour reversal during
the minimums in brightness.

We can confirm our conclusions, made in our previous paper (Semkov et al.
[7]), that the photometric properties of GM Cep can be explained by a
superposition of highly variable accretion from the circumstellar disk onto
the stellar surfice and occultation from circumstellar clumps of dust,
planetesimals or from other features of the circumstellar disk.

## Acknowledgements

This work was partly supported by the Bulgarian Scientific Research Fund of
the Ministry of Education and Science under the grants DN 08-1/2016, DN
18-10/2017 and DN 18-13/2017.

## References

  * [1] Chen, W. P., Hu, S. C.-L., Errmann, R., et al., 2012, ApJ 751, 118
  * [2] Chen, W. P., and Hu, S. C.-L., 2014, IAUS 293, 74
  * [3] Contreras, M. E., Sicilia-Aguilar, A., Muzerolle, J., et al. 2002, AJ, 124, 1585
  * [4] Huang P. C., Chen W. P., Mugrauer M., Bischoff R., et al., 2019, ApJ 871, 183
  * [5] Marschall, L. A. and van Altena, W. F., 1987, AJ, 94, 71
  * [6] Semkov, E. H., and Peneva, S. P., 2012, ApSS 338, 95
  * [7] Semkov, E. H., Ibryamov, S. I., Peneva, S. P., et al., 2015, PASA 32, e011
  * [8] Sicilia-Aguilar, A., Hartmann, L., Hernandez, J., Briceno, C., and Calvet, N., 2005, AJ 130, 188
  * [9] Sicilia-Aguilar, A., Merın, B., Hormuth, F. et al., 2008, ApJ 673, 382
  * [10] Xiao, L., Kroll, P., and Henden, A., 2010, AJ 139, 1527