# Community Detection From Multiple Observations: From Product Graph Model to Brain Applications

Tiziana Cattai, Gaetano Scarano, Marie-Constance Corsi, Fabrizio Devico

Fallani, Stefania Colonnese

<figure><img src="image_1.png" /></figure>

## To cite this version:

Tiziana Cattai, Gaetano Scarano, Marie-Constance Corsi, Fabrizio Devico Fallani, Stefania Colonnese. Community Detection From Multiple Observations: From Product Graph Model to Brain Applications. IEEE Transactions on Signal and Information Processing over Networks, 2025, 11, pp.201-214. 10.1109/TSIPN.2025.3540702. hal-05273024

HAL Id: hal-05273024

https://inria.hal.science/hal-05273024v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Community Detection from Multiple Observations: from Product Graph Model to Brain Applications

Tiziana Cattai, Member, IEEE, Gaetano Scarano, Marie-Constance Corsi, Fabrizio De Vico Fallani, Stefania Colonnese, Senior Member, IEEE

***Abstract***—This paper proposes a multilayer graph model for community detection based on multiple observations. This scenario is common when different estimators are used to infer graph edges from signals at the nodes, or when various signal measurements are taken. The multilayer network stacks these graph observations at different layers and links replica nodes at adjacent layers. This configuration corresponds to the Cartesian product between the ground truth graph and a path graph, where the number of nodes matches the number of observations. Using the algebraic structure of the Laplacian of the Cartesian multi-layer network, we infer a subset of the eigenvectors of the true graph and perform community detection. Experimental results on synthetic graphs demonstrate the accuracy of the method, which outperforms state-of-the-art approaches in correctly detecting graph communities. Finally, we apply our method to distinguish between different brain networks derived from real EEG data collected during motor imagery experiments. We conclude that our approach is promising for identifying graph communities when multiple graph observations are available, and it shows potential for applications such as EEG-based motor imagery applications.

***Index Terms***—Multilayer Networks, Product Graph, Community Detection, Functional Connectivity.

## I. INTRODUCTION

Various cutting-edge fields, such as computational biology, precision medicine, and social studies, collect measurements over the nodes of a network, thus driving the development of new tools of graph signal processing (GSP) [1], [2]. Extracting the community structure of a network from these observations enhances the analysis and understanding of the data [3], [4].

The human brain stands out as a perfect challenge for GSP-based studies. In this context, researchers represent brain data collected at different EEG electrodes as Signals on Graphs (SoG) [1], [2], modeling the brain's structure with graphs. GSP methods can identify network communities and facilitate brain network representation [5] or mental state tracking [6] highlighting interactions between different brain regions through the use of different recording techniques. While EEG signals are often used for their practicality, ease of measurement, and low cost, other types of signals, such as functional magnetic resonance imaging (fMRI), have also been used in the field of graph signal processing. GSP methods have been effectively applied to fMRI data to analyze the dependencies between

brain region activations and their physical connections [5]. Moreover, graph sampling techniques in GSP provide an alternative to traditional component reduction approaches for classification tasks [7], and properties of the Graph Fourier Transform from fMRI time series have proven useful in such contexts [8].

EEG-based studies infer brain interactions using various connectivity estimators between nodes [9] or the same estimator applied across different frequency bands, [10], each yielding distinct brain network representations [11].

In this and many other applications, various graph learning methods produce multiple estimates of a single unknown ground truth network. This raises an important question: can we combine these multiple observations to detect the true communities with higher accuracy?

In this work, i) we introduce a novel model for handling multiple estimates of a single unknown ground truth graph by considering them as different layers of a multilayer product graph [12], which depend on both the unknown ground truth graph and the number of estimates [13]; ii) we leverage this multilayer model to estimate the ground truth eigenvectors [14], [15] and detect network communities by appropriately scaling and clustering these eigenvectors.

Several state-of-the-art methods identify communities on multilayer graphs through graph structure analysis [16], [17], graph matrix factorization [18], [19], pattern mining functions [20], or specific assumptions on the signal on graph [21]. Unlike these methods [6], we leverage the algebraic structure of the multilayer graph to model the redundancy in multiple observations of the underlying ground truth graph. This approach enables us to fully exploit the multiple observations to accurately detect the true communities. Results on synthetic data show that the proposed method outperforms state-of-the-art methods [22], [23], [21], [24], [25]. Furthermore, we present experimental results on real brain data and show that integrating different graph estimates into the multilayer model significantly improves community detection performances.

To sum up, the paper addresses the challenge of community detection from different estimates of an unknown underlying graph, with a particular focus on brain data analysis. The main contributions are as follows:

* we propose modeling multiple, possibly noisy observations of the ground truth graph as a multilayer network;

* we adopt a Cartesian Product Graph Model (CPGM) for the multilayer graph, constructed as the Cartesian product graph of the underlying ground truth graph and

T. Cattai, S. Colonnese, G. Scarano are with Dept. of Information Engineering, Electronics and Telecommunication, Sapienza University of Rome, Italy.

F. De Vico Fallani and M.C. Corsi are with Sorbonne Université, Paris Brain Institute - ICM, CNRS, Inria, Inserm, AP-HP, Hopital Pitie Salpetriere, Paris, France.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Graphical abstract of the proposed community detection framework. From multiple adjacency matrices $A^{(m)}$, the CPGM generates a combined adjacency matrix $\mathcal{A}$ and the associated Laplacian matrix $\mathcal{L}$. Then, eigenvectors $\mathcal{U}$ are computed and the low-frequency components are firstly rearranged and flipped and then used for spectral clustering.</figcaption></figure>

a path graph; we leverage the algebraic CPGM structure to estimate the eigenvectors of the ground truth graph;

* we present a novel CPGM-based Community Detection (CPGM-CD) algorithm of the ground truth graph, which appropriately scales the estimated first eigenvectors of the ground truth graph and applies spectral clustering;

* after demonstrating that CPGM-CD outperforms various state-of-the-art algorithms, we apply it to real EEG data recorded during motor imagery-based Brain Computer Interface (BCI) experiments. The CPGM, built on several state-of-the-art functional connectivity estimators, utilizes the diversity in deriving the brain network and captures different kinds of interactions. We show that CPGM-CD effectively highlights relevant communities, paving the way for a deeper understanding of human brain behaviors.

In Fig. 1 we present an overview of our method.

The structure of the paper is as follows. Section II describes the observation model. In Section III we introduce the Cartesian Product Graph Model (CPGM) used to merge multiple graph observations. Section IV presents the CPGM-based community detection (CPGM-CD) algorithm. In Section V, we report results on synthetic data, after the description of the generation data model. In Section VI we apply the proposed method to real brain EEG signals, and Section VII concludes the paper.

In Table I, we list the main notation used in the paper.

## II. OBSERVATION MODEL

The unknown ground truth graph $\mathcal{G}_{GT}=(\mathcal{V}_{GT}, \mathcal{E}_{GT})$ is defined as a set of $N$ nodes (or vertices) and links (or edges). The ground truth adjacency matrix $\mathbf{A}_{GT} \in \mathbb{R}^{N \times N}$ element $a_{ij}^{(GT)}$ represents the weight of the link connecting the nodes $i$ and $j$. If the ground truth graph $\mathcal{G}_{GT}$ is endowed with a community structure, the vertex set $\mathcal{V}_{GT}$ is partitioned into

$P$ disjoint sets, referred to as node communities. The $n$-th vertex has an associated community label $\gamma_n \in \{0, ..., P-1\}$. The community structure is represented by a $P \times N$ affinity matrix $\Gamma$ whose $n$-th column is a unit-norm binary vector $1_p \in \{0, 1\}^P$ whose $p$-th component equals to 1 if the $n$-th node belongs to the $p$-th community and zero otherwise.

The graph communities may reflect different network structures [26], [27] based on the conditioning of the affinity matrix $\Gamma$ on the adjacency matrix $\mathbf{A}$. Herein, we restrict ourselves to the case of assortative networks, which is by far the most common in the applications, where nodes belonging to the same community are more likely to be linked than nodes belonging to different communities [28]. Let us remark that widely adopted random network models, such as the Stochastic Block Models (SBM) and the SBM subset of Planted Partition Models (PPM), are assortative networks [29], [30].

The graph Laplacian $\mathbf{L}_{GT} \in \mathbb{R}^{N \times N}$ is defined as $\mathbf{L}_{GT} = \mathbf{D}_{GT} - \mathbf{A}_{GT}$ where $\mathbf{D}_{GT}$ is the diagonal degree matrix, with diagonal entries $d_i = \sum_{j=1}^{N} a_{ij}^{(GT)}$, $i=0, ..., N-1$. Important features in Graph Signal Processing (GSP) theory are associated to the eigendecomposition of the graph Laplacian $\mathbf{L}_{GT}$:

$$ \mathbf{L}_{GT} = \mathbf{U}_{GT} \mathbf{\Lambda}_{GT} \mathbf{U}_{GT}^{T} \quad (1) $$

where the diagonal matrix $\mathbf{\Lambda}_{GT}$ and the orthonormal matrix $\mathbf{U}_{GT}$ contain the Laplacian eigenvalues $\lambda_n^{(GT)}$ and eigenvectors $\mathbf{u}_n^{(GT)}$, $n=0, ..., N-1$, respectively. The $n$-th Laplacian eigenvalue $\lambda_n^{(GT)}$ reflects the properties of the $n$-th eigenvector $\mathbf{u}_n^{(GT)}$ [31]: eigenvectors with large eigenvalues are more robust to noise [32], [33]; eigenvectors with small eigenvalues change smoothly on connected nodes, and carry information on graph connectivity [34], [35].

The unknown ground truth graph adjacency matrix $\mathbf{A}_{GT}$ is accessible via $M$ observations, modeled as follows:

$$ \mathbf{A}^{(m)} = \mathbf{A}_{GT} + \mathbf{W}^{(m)}, \quad m = 0, \dots, M-1 \quad (2) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Cartesian product graph model of multiple observations: example of $M = 3$ observations of a ground truth graph with $N = 5$ nodes (first row), of the multilayer network built on them (second row), and of the Cartesian product graph model (CPGM) of the observations. The CPGM is obtained from the ground truth graph $\mathcal{G}_{GT}$ and the $M$ nodes path graph $\mathcal{G}_M$.</figcaption></figure>

where $\mathbf{A}_{GT}$ is a random matrix drawn from a conditional distribution

$$ \mathbf{A}_{GT} \approx p_{\mathbf{A}|\Gamma}(\mathbf{A}|\Gamma), \quad (3) $$

and $\mathbf{W}^{(m)}$, $m = 0, \dots, M-1$ are i.i.d. random symmetric matrices independent on $\mathbf{A}_{GT}$ modeling the fluctuations -due to sample noise or estimation errors- of the observed matrices $\mathbf{A}^{(m)}$ with respect to the ground-truth one $\mathbf{A}_{GT}$.

We leverage this model to develop a community detection algorithm from multiple observations.

### III. FROM MULTIPLE OBSERVATIONS TO CARTESIAN PRODUCT GRAPH MODEL (CPGM)

State-of-the-art methods may learn the graph from node signals with different criteria, resulting in various observed adjacency matrices [36], [37].

Herein, we obtain the $M$ observations of the graph using a multilayer network model [12] defined as follows:

* the multilayer network consists of $M$ layers, each containing $N$ nodes;

<table><caption>TABLE I
TABLE OF MAIN NOTATION.</caption><thead><tr><th>Notation</th><th>Description</th></tr></thead><tbody><tr><td>$\mathcal{A}, \mathbf{A}_{GT}, \mathbf{A}_M$</td><td>adjacency matrix of multilayer, ground truth and path graph</td></tr><tr><td>$\mathcal{V}, \mathcal{N}$</td><td>set and number of nodes at each layer</td></tr><tr><td>$\mathcal{E}, N_c$</td><td>set and number of links at each layer</td></tr><tr><td>$\mathcal{L}, \mathbf{L}_{GT}, \mathbf{L}_M$</td><td>Laplacian matrix of multilayer, ground truth and path graph</td></tr><tr><td>$\mathcal{U}, \mathbf{U}_{GT}, \mathbf{U}_M$</td><td>Eigenvectors of multilayer network, ground truth and path graph</td></tr><tr><td>$\mathbf{M}, \mathbf{M}-1$</td><td>number of layers and links between corresponding nodes in different layers</td></tr><tr><td>$N_{CD}$</td><td>number of selected eigenvectors</td></tr><tr><td>$\tilde{\mathbf{U}}, \tilde{\mathbf{U}}$</td><td>matrix of selected eigenvectors, normalized matrix of selected eigenvectors</td></tr><tr><td>$\mathbf{P}$</td><td>number of communities</td></tr></tbody></table>

* the estimated adjacency matrix for the $m$-th layer equals to the $m$-th estimated adjacency matrix, where $m = 0, \dots, M-1$;

* the corresponding nodes between consecutive layers' pairs $(m, m+1)$, $m = 0, \dots, M-2$ are connected.

The adjacency matrix $\mathcal{A}$ of the multilayer network is therefore constructed with $MN \times N$ blocks on the diagonal, representing the intra-layer links, along with $M-1$ upper-diagonal and $M-1$ lower-diagonal blocks equal to the $N \times N$ identity matrix. Fig. 2 (first row) shows the multilayer network model for the case of $M=3$ observations of a ground truth graph of $N=5$ nodes.

Let us remark that for ideal observations, the multilayer network boils down to $M$ replicas of the ground truth graph, with links between corresponding nodes.

Our proposed model corresponds to a specific kind of multilayer graph [38], where the network is the Cartesian product graph of the ground truth graph $\mathcal{G}_{GT}$ and a length $M$ path graph $\mathcal{G}_M$:

$$ \mathcal{G} = \mathcal{G}_{GT} \square \mathcal{G}_M \quad (4) $$

The Cartesian product graph model has been applied [38] under the assumption that the observed graphs are perturbed versions of the same graph, and that the signal at each node is correlated across different layers. In our case, this reflects the assumption that different estimators observe the same underlying functional connectivity, resulting in correlations between them. Thus, we model the $M$ observations as layers of a Cartesian product graph between $\mathcal{G}_{GT}$ and $\mathcal{G}_M$. Fig. 2 (second row) shows an example of the Cartesian Product Graph Model (CPGM) of multiple observations, where $M = 3$ observations of a ground truth graph with $N = 5$ nodes are available (first row). The multilayer network model (second row) corresponds to the Cartesian product graph between the ground truth graph $\mathcal{G}_{GT}$ and the $M$ nodes path graph $\mathcal{G}_M$ (third row).

The Cartesian Product Graph Model (CPGM) endows the observations with a strong algebraic structure. Indeed, the adjacency matrix of the Cartesian product graph $\mathcal{G} = \mathcal{G}_{GT} \square \mathcal{G}_M$, results as the Kronecker sum of the adjacency matrices of $\mathbf{A}_{GT}$ and $\mathbf{A}_M$:

$$ \mathcal{A} = \mathbf{A}_{GT} \oplus \mathbf{A}_M \quad (5) $$

with $\mathbf{A}_M$ as the adjacency matrix of a path graph of order $M$. The graph Laplacian is obtained as the following Kronecker product:

$$ \mathcal{L} = \mathbf{L}_{GT} \otimes \mathbf{L}_M \quad (6) $$

Given the eigenvectors matrices $\mathbf{U}_{GT}, \mathbf{U}_M$ of the ground truth $\mathcal{G}_{GT}$ and path graph $\mathcal{G}_M$, respectively, the CPGM eigenvectors are computed as follows [38]:

$$ \mathcal{U} = \mathbf{U}_M \otimes \mathbf{U}_{GT} \quad (7) $$

Besides, from the literature in [38] we recognize that the CPGM eigenvalues $\lambda_k, k = 0, \dots, M \cdot N - 1$ are computed by
<figure><img src="image_5.png" /><figcaption>Fig. 3. Example of eigenvalues' pairs ($\lambda_m^{(M)}$, $\lambda_n^{(GT)}$), for $M = 3$, $N = 6$; the corresponding CPGM eigenvalues $\lambda_k$ are obtained by summing the components $\lambda_m^{(M)} + \lambda_n^{(GT)}$. Let us observe that the order of the eigenvalues is independent from the ordering of the layer, and it reflects the order of the eigenvalues sum. Geometrically, the value of the eigenvalues' sum is indicated by the orange slanting lines, where the arrow indicates increasing sum levels. All pairs $\lambda_m^{(M)}$, $\lambda_n^{(GT)}$ that lie below a certain slanted line correspond to CPGM eigenvalues $\lambda_k$ that are below a given threshold. The ($\lambda_m^{(M)}$, $\lambda_n^{(GT)}$) pairs in the orange area contribute to the smallest $N_{CD}$ CPGM eigenvalues, and involve a few of the smallest $N_{CD}$ eigenvalues $\lambda_{n_k}^{(GT)}$, $k = 0, ..., N_{CD}$ of the unknown ground truth Laplacian $L_{GT}$.</figcaption></figure>

summing up each and every pair of eigenvalues $\lambda_n^{(GT)}$, $n = 0, \dots, N-1$, and $\lambda_m^{(M)}$, $m=0, \dots, M-1$, as follows:

$$ \lambda_k \in \left\{ \begin{array}{l} \lambda_m^{(M)} + \lambda_n^{(GT)}, \\ m=0, \dots, M-1, n=0, \dots, N-1 \end{array} \right\} \quad (8) $$

This is illustrated in Fig.3, where the x- and y- axes represent $\lambda_n^{(GT)}$ and $\lambda_m^{(M)}$, respectively. Each point is associated to a multilayer eigenvalue, associated to the sum ($\lambda_m^{(M)} + \lambda_n^{(GT)}$), for $M = 3$, $N = 6$. Let us observe that the points below the orange sloped curves with x-axis intersect $\bar{\lambda}$ are associated to multilayer eigenvalues smaller than $\bar{\lambda}$.

The $NM \times NM$ CPGM eigenvector matrix $\mathcal{U}$ obtained by the Kronecker product of $U_M$ and $U_{GT}$ is endowed with a strong structure: it stores $M$ differently scaled replicas of the ground truth eigenvectors $\mathbf{u}_n^{(GT)}$, with $n = 0, \dots, N-1$, with scaling factors equal to the coefficients of the eigenvectors in $U_M$.

Both the eigenvalues $\lambda_m^{(M)}$ and eigenvectors $U_M$ of the length $M$ path graph $\mathcal{G}^{(M)}$ are known in the literature, namely

$$ \lambda_m^{(M)} = 2 \left( 1 - \cos \left( \frac{\pi m}{M} \right) \right) \quad (9) $$

and

$$ \mathbf{u}_m^{(M)}[l] = \cos \left( 2\pi \frac{(m-1/2)}{2M} l \right). \quad (10) $$

where $m=0, \dots, M-1$ is the path graph eigenvalue index, and the variable $l=0, 1, \dots, M-1$ denotes the $l$-th coefficient of the $m$-th path graph eigenvector.

Fig. 4 shows a toy example ($M = 3$, $N = 4$) of the Cartesian product graph $\mathcal{G}$ between the ground truth $\mathcal{G}_{GT}$ and the path graph $\mathcal{G}_M$, of adjacency and eigenvectors matrices, $\mathbf{A}_{GT}$, $\mathbf{A}_M$, and $\mathbf{U}_{GT}$, $\mathbf{U}_M$, respectively. Fig. 4 also depicts the CPGM adjacency and the eigenvectors matrices $\mathcal{A}$, $\mathcal{U}$. The CPGM eigenvector matrix $\mathcal{U}$ contains $M^2$ replicas of the $n$-th

<figure><img src="image_6.png" /><figcaption>Fig. 4. CPGM algebraic structure for a graph $\mathcal{G} = \mathcal{G}_{GT} \square \mathcal{G}_M$ with $M = 3$ layers and $N = 4$ nodes each. The CPGM adjacency matrix $\mathcal{A}$ is given by $\mathbf{A}_{GT} \oplus \mathbf{A}_M$. The CPGM eigenvector matrix $\mathcal{U}$ is the Kronecker product $\mathbf{U}_{GT} \otimes \mathbf{U}_M$ of the path and ground truth graphs' eigenvectors $\mathbf{U}_{GT}$ and $\mathbf{U}_M$. Therefore, $\mathcal{U}$ contains $M^2$ replicas of the $n$-th ground truth eigenvector $\mathbf{u}_{GT}^{(n)}$, $n = 0, \dots, N-1$, suitably displaced and scaled. The scaling factor and position of each replica of $\mathbf{u}_{GT}^{(n)}$ depend on the ground truth eigenvalues and are unknown.</figcaption></figure>

ground truth eigenvector $\mathbf{u}_n^{(GT)}$, $n = 0, \dots, N-1$, displaced and scaled.

Let us remark that the positions of the replicas of the $n$-th eigenvector $\mathbf{u}_{GT}^{(n)}$ in $\mathcal{U}$ depend on the ordering of the eigenvalues $\{\lambda_n^{(GT)} + \lambda_m^{(M)}, m=0, \dots, M-1, n=0, \dots, N-1\}$, i.e., on the unknown eigenvalues $\lambda_n^{(GT)}$. This prevents direct estimation of the eigenvectors $\mathbf{u}_{GT}^{(n)}$ from the multilayer network eigenvectors matrix $\mathcal{U}$.

In the following Section, we propose a community detection (CD) algorithm that leverages the partial knowledge on the CPGM eigenvector matrix $\mathcal{U}$ to detect ground truth communities.

## IV. CPGM-BASED COMMUNITY DETECTION (CPGM-CD)

The CPGM eigenvectors matrix $\mathcal{U}$ gathers different replicas of the ground truth eigenvectors $\mathbf{u}_{GT}^{(n)}$, $n = 0, \dots, N-1$. We propose a novel CPGM-based Community Detection (CPGM-CD) methodology that exploits this algebraic structure by

1. extraction of $N_{CD} < N$ eigenvectors from the CPGM eigenvector matrix $\mathcal{U}$;

2. scaling the replicas to obtain a set $\tilde{\mathcal{U}}_{CD}$ of estimated ground truth eigenvectors;

3. clustering on the estimated eigenvectors' set $\mathcal{U}_{CD}$, assigning a community label $\gamma_n$ to each node.

### A. Extraction Stage

The first stage extracts the first $N_{CD}$ eigenvectors in $\mathcal{U}$, ordered for ascendant eigenvalues $\lambda_k$, $k=0, \dots, N_{CD} - 1$. The following $NM \times N_{CD}$ matrix is built:

$$ \mathcal{U}_{CD} = [\mathbf{u}_0 \cdots \mathbf{u}_{N_{CD}-1}] \quad (11) $$

The first $N_{CD}$ eigenvectors represent the low frequency components of the Graph Fourier Transform (GFT) basis for the multilayer network, and tend to smoothly vary over connected areas [1].
Each eigenvector $\mathbf{u}_k, k = 0, \dots, N_{CD}$ contains $M$ replicas of a ground truth eigenvector $\mathbf{u}_{n_k}^{(GT)}$. To determine whether the ground truth eigenvectors $\mathbf{u}_{n_k}^{(GT)}, k = 0, \dots, N_{CD}$ are low frequency or not, we observe that the CPGM eigenvalues $\lambda_k, k = 0, \dots, MN-1$ are computed by summing up each and every pair of eigenvalues $\lambda_n^{(GT)}, n = 0, \dots, N-1$, $\lambda_m^{(M)}, m=0, \dots, M-1$. Hence, the $N_{CD}$ smallest eigenvalues $\lambda_k, k=0, \dots, N_{CD}-1$ correspond to the $N_{CD}$ pairs $(m, n)$ with the smallest sums $\lambda_n^{(GT)} + \lambda_m^{(M)}$.

Fig.3 shows an example of pairs $(\lambda_m^{(M)}, \lambda_n^{(GT)})$, for $M = 3, N = 6$; the corresponding eigenvalues $\lambda_k, k=0, \dots, MN-1$ are ordered in terms of increasing sum. The extraction stage selects the eigenvectors corresponding to the smallest $N_{CD}$ CPGM eigenvalues, represented by the $(\lambda_m^{(M)}, \lambda_n^{(GT)})$ pairs in the orange area. The indexes $n_k, k = 0, \dots, N_{CD}-1$ of the ground truth eigenvalues $\lambda_{n_k}^{(GT)}, k = 0, \dots, N_{CD}$ contributing to the smallest $N_{CD}$ eigenvalues $\lambda_k, k = 0, \dots, N_{CD}-1$ depend on the unknown ground truth Laplacian $L_{GT}$. Still, it is clear that the set of $\lambda_{n_k}^{(GT)}, k = 0, \dots, N_{CD}$ include some of the smallest eigenvalues $\lambda_n^{(GT)}$. Hence, the ground truth eigenvectors $\mathbf{u}_{n_k}^{(GT)}, k = 0, \dots, N_{CD}-1$ appearing in the first $N_{CD}$ eigenvectors of $\mathcal{U}$ represent low frequency GFT components of the ground truth graph.

As shown in Fig.4, the matrix $\mathcal{U}_{CD} \in \mathbb{R}^{MN \times N_{CD}}$ of the first $N_{CD}$ eigenvectors $\mathbf{u}_k, k=0, \dots, N_{CD}-1$ collects $M \times N_{CD}$ scaled replicas

$$ \mathbf{s}_m^{(k)} = \alpha_m^{(n_k)} \mathbf{u}_{n_k}^{(GT)} \quad (12) $$

of the low frequency ground truth eigenvectors $\mathbf{u}_{n_k}^{(GT)}, k = 0, \dots, N_{CD}$, where the unknown factor $\alpha_m^{(n_k)}$ accounts for the multiplying coefficient in (7).

We denote by $S_{CD}$ the $N \times MN_{CD}$ matrix defined as

$$ S_{CD} = [\mathbf{s}_0^{(0)} \cdots \mathbf{s}_{M-1}^{(N_{CD}-1)}]. \quad (13) $$

Fig.5 depicts, from left to right, the matrices $\mathcal{U} \in \mathbb{R}^{MN \times MN}$, $\mathcal{U}_{CD} \in \mathbb{R}^{MN \times N_{CD}}$ and $S_{CD} \in \mathbb{R}^{N \times MN_{CD}}$ for the toy case of a multilayer network with $M = 3$ layers and $N = 6$ nodes at each layer.

We emphasize that the choice of the parameter $N_{CD}$ determines the number of eigenvectors used for the clustering. In this study, we build on results obtained in [31] and select the 20% of eigenvectors associated to the smallest eigenvalues. For proper clustering, the column vectors in $S_{CD}$ are scaled as described below.

## B. Scaling Stage

The columns of the matrix $S$ are proportional, but not identical, to the low-frequency eigenvectors of the ground truth graph. We apply a scaling procedure that: i) normalizes the vectors' energy, and ii) flips the vectors so as to align the sign of their first coefficient.

This is achieved by computing the normalized matrix for clustering as follows:

$$ \tilde{\mathcal{U}}_{CD} = \mathbf{E}^{-1/2} \cdot \mathbf{B} \cdot S_{CD} \quad (14) $$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Graphical representation of eigenvectors of the multilayer network with $M = 3$ layers and $N = 6$ nodes at each layer. In panel (a) we report the eigenvector matrix $\mathcal{U}$. In (b) we show that we select the first $N_{CD}$ eigenvectors obtaining $\mathcal{U}_{CD}$ with $N \times M = 18$ rows and $N_{CD}$ columns. Then, we extract the vectors and arrange them into the matrix $\tilde{\mathcal{U}}_{CD}$ of dimensions $N \times N_{CD}M$ as in panel (c). Note that in this specific configuration, there is no need of sign-flipping for the selected eigenvectors and $\tilde{\mathcal{U}}_{CD} = S_{CD}$.</figcaption></figure>

with $\mathbf{E}$ as the diagonal energy matrix:

$$ \mathbf{E} = \begin{bmatrix} \|\mathbf{s}_0^{(0)}\|_2 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & \|\mathbf{s}_{M-1}^{(N_{CD}-1)}\|_2 \end{bmatrix} \quad (15) $$

and $\mathbf{B}$ the diagonal binary flipping matrix

$$ \mathbf{B} = \begin{bmatrix} \text{sign}(\mathbf{s}_0^{(0)}[0]) & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & \text{sign}(\mathbf{s}_{M-1}^{(N_{CD}-1)}[0]) \end{bmatrix} \quad (16) $$

For clarity, we emphasize that the proposed scaling stage operates independently of the clustering. Although the joint design of the scaling and clustering steps could be promising, it is out beyond the scope of this paper and is left for future work.

## C. Clustering Stage

Community detection can be achieved by clustering of the Laplacian eigenvectors, a method known as spectral clustering [39], [40]. In classical spectral clustering, the $i$-th node of the ground truth graph is associated with the $i$-th row of the eigenvector matrix $U_{GT}$. Clustering is then performed on the rows of $U_{GT}$, with nodes whose rows fall within the same cluster assigned the same community label. Herein, we perform clustering over the $N$ rows of $\tilde{\mathcal{U}}_{CD}$, which contains $N_{CD}$ scaled replicas of the ground truth low frequency eigenvectors, and assign a community label $\gamma_n$ to each node.

To sum up, recent literature on signal-on-graphs has shown that eigenvectors associated to small eigenvalues typically describe the community structure of a graph [1], [2]. Following the same strategy, we extend this concept to the multilayer network by using the eigenvectors associated with the smallest
eigenvalues of the multilayer graph for community detection. After computing the first $N_{CD}$ eigenvectors of $\mathcal{U}$, we rearrange them into the matrix $S_{CD}$ of size $N \times MN_{CD}$. We then regularize the energy and direction of each eigenvector estimate by scaling each column of $S_{CD}$ to unitary energy and flipping the sign so that its first element is positive. Finally, we identify the communities through spectral clustering of the regularized eigenvectors in $\tilde{\mathcal{U}}_{CD}$ [41], [42], [43]. The spectral clustering provides a community label for each node of the ground truth graph.

Remarkably, CPGM-CD acts on the multilayer network of $M$ layers representing multiple observations and identifies the communities of the single underlying ground truth graph. The sequence of main steps of CPGM-CD is presented in Algorithm 1.

Let us remark that one of the known parameters of our method is the number of graph communities. Once we have estimated the ground-truth eigenvectors $U_{GT}$, we perform spectral clustering with a selected number of clusters to be identified. The number of communities can be estimated leveraging state-of-the-art clustering techniques. Many studies have explored ways to determine the optimal number of clusters for a given data structure, such as the well established consensus clustering [44], or by using methods that automatically find the optimal number of clusters, such as silhouette analysis [45] or the Calinski-Harabasz index [46], [47]. Specifically, computing the consensus matrices that indicate the persistence of clustering across different layers can help estimate the number of communities. For example, this can be achieved by selecting the number of communities that yields the most stable or persistent, communities across different layers. The consensus matrix across layers can also be studied as a feature for graph classification [48]. Alternatively, the eigengap method can be used to estimate the number of communities [49]. An in-depth analysis of these topics is beyond the scope of this paper and is left for future work.

## V. EXPERIMENTAL RESULTS

Herein, we evaluate the CPGM-CD method on stochastic graphs, showing that it outperforms state-of-the-art competitors, particularly in cases where the observations are affected by large fluctuations.

We generate the $M$ observations according to a Stochastic Block Model, known as the Planted Partition Model (PPM), where the set of observable links is related to the unknown community structure $\Gamma$ as follows. The $i, j$ element $A_{ij}$ of the ground truth matrix $\mathbf{A}_{GT}$ is drawn from the following conditional distribution:

$$ p_{A_{ij}|\gamma_i, \gamma_j}(A_{ij}|\gamma_i, \gamma_j) = \begin{cases} p_s\delta(A_{ij}) + (1-p_s)\delta(A_{ij}-1) & \gamma_i = \gamma_j \\ p_d\delta(A_{ij}) + (1-p_d)\delta(A_{ij}-1) & \gamma_i \neq \gamma_j \end{cases} \quad (17) $$

where the probabilities of observing a link between nodes of the same community or across different communities are denoted as $p_s$ and $p_d$, respectively.

To evaluate the method's performance in the presence of connectivity estimation errors, we introduce additive Gaussian noise. This approach is justified by the characteristics of brain

## Algorithm 1 CPGM-CD Algorithm

Input: $\mathcal{A}, \tilde{N}_{CD}$

Output: Community partition $\Gamma$

1: Step 1: Eigendecomposition of $\mathcal{L}$

a: Compute the Graph Laplacian of the multilayer network

$$ \mathcal{L} = \mathcal{D} - \mathcal{A} $$

b: Compute the eigendecomposition

$$ \mathcal{L} = \sum_{q=0}^{N \times M - 1} \lambda_q \mathbf{u}_q \mathbf{u}_q^H $$

2: Step 2: Community Detection

a: Select the first $N_{CD}$ eigenvectors of $\mathcal{L}$: $\tilde{\mathcal{U}}_{CD} = [\mathbf{u}_0 \cdots \mathbf{u}_{N_{CD}-1}]$

b: Reorganize $\tilde{\mathcal{U}}$ in a strip of $N$ rows and $N_{CD} \times M$ columns

c: Normalize the eigenvectors by their energy and force them to be in the same direction:

$$ \tilde{\mathcal{U}}_{CD} = \mathbf{E}^{-1/2} \cdot \mathbf{B} \cdot S_{CD} $$

d: $\Gamma \leftarrow$ Spectral clustering $(\tilde{\mathcal{U}}_{CD})$

signals, where connectivity estimators are typically based on the cross-correlation of signals recorded over time at the electrodes [50], [51]. Since multiple components are combined, the fluctuations of the estimator are well approximated by a bell-shaped (Gaussian) distribution [52]. Therefore, to validate the algorithm in a controlled synthetic scenario, we consider additive Gaussian noise. Specifically, we consider the case where the $M$ observations differ due to normally distributed estimation error, i.e.,

$$ \mathbf{A}^{(m)} = \mathbf{A}_{GT} + \mathbf{W}^{(m)}, \quad m = 0, \dots, M-1 \quad (18) $$

where $\mathbf{W}^{(m)}, m = 0, \dots, M-1$ are i.i.d. Gaussian matrices at Signal-to-Noise-Ratio SNR¹.

Fig. 6 refers to the toy case of a ground truth graph of $N=4$ nodes belonging to $P=2$ communities. We considered $M=3$ observations, generated according to the PPM model using $p_s = 0.9, p_d = 0.1$ and the ground truth graph $\mathcal{G}_{GT}$ is sketched in Fig. 6 (a), where the node color represents the node ground truth community label. The $M=3$ observations $\mathbf{A}^{(m)}, m = 0, 1, 2$ are affected by an additive Gaussian noise at SNR=40 dB. Then, a multilayer network formed by $M=3$ layers with $N=4$ nodes each is built by collecting $M, N \times N$ adjacency matrices. Specifically, Fig. 6 (b) plots the related multilayer adjacency matrix $\mathcal{A}$, where the diagonal blocks represent the observations $\mathbf{A}^{(m)}, m = 0, 1, 2$, and the upper and lower diagonal blocks are identity matrices representing the links between corresponding nodes on adjacent layers.

$$ w_{ij}^{(m)} = w_{ji}^{(m)} \approx \mathcal{N}(0, \sigma_w^2/2), \text{ where} $$

$$ \sigma_w^2 = \frac{1}{N^2} \|\mathbf{A}_{GT}\|_F^2 \cdot \text{SNR}^{-1} $$

and $\| \cdot \|_F$ denotes the matrix Frobenius norm.
<figure><img src="image_8.png" /><figcaption>Fig. 6. (a) Graph with $N = 4$ nodes belonging to $P = 2$ communities generated according to (17), (18) with $p_s = 1 - p_d = 0.9$ and SNR = 40 db; the node color represents its community label; (b) multilayer adjacency matrix $\mathcal{A}$ built on $M = 3$ observations $\mathbf{A}^{(m)}$; (c) ground truth and path graph eigenvector matrices $U_{GT}$, $U_M$, and ground truth and estimated multilayer eigenvector matrices $\mathcal{U} = U_{GT} \otimes U_P$, $\hat{\mathcal{U}}$. The CPGM-CD algorithm exactly recovers the community partition represented by different colors in panel (a).</figcaption></figure>

The algebraic structure of the multilayer adjacency matrix $\mathcal{A}$ reflects into the structure of the multilayer eigenvector matrix $\mathcal{U}$.

Fig. 6 (c) represents the $4 \times 4$ ground truth eigenvector matrices $U_{GT}$, the $3 \times 3$ path graph eigenvector matrices $U_M$, the $12 \times 12$ multilayer ground truth eigenvector matrices $\mathcal{U}$ -obtained as the Kronecker product of $U_{GT}$ and $U_M$, and the estimated multilayer eigenvector matrix $\hat{\mathcal{U}}$ -obtained by eigenanalysis of the multilayer Laplacian $\mathcal{L}$ obtained from $\mathcal{A}$. The coefficients of the first estimated eigenvectors in $\hat{\mathcal{U}}$ in Fig. 6 (c), apart from the different ordering and from the energy normalization which reflect into different colors in the plot, match the first eigenvectors based on the CP model in $\mathcal{U}$. These vectors are selected and reshaped to apply the scaling and clustering stages of the CPGM-CD algorithm that correctly recovers the nodes community labels $\gamma_i$, $i = 0, \dots, 3$.

Fig. 7 plots the result of the community detection algorithm on three PPM graphs of $N = 45$ nodes, belonging to $P = 2, 3, 4$ communities, and generated according to $p_s=0.8, p_d=0.2$. The $M = 3$ observations are obtained at SNR = 40. The CPGM-CD algorithm is applied by leveraging $N_{CD} = 3$ eigenvectors, of $NM = 45 \cdot 3$ coefficients each, suitably reshaped into the $N \times M \cdot N_{CD} = 45 \times 9$ matrix $S_{CD}$. After proper scaling and clustering, the CPGM-CD method correctly recovers the community partition, identified by different colors in Fig.7.

We now compare the performance of the CPGM-CD method with those achieved by: i) separately applying state-of-the-art community detection algorithms to $\mathbf{A}^{(m)}, m=0, \dots, M-1$ and ii) assigning the community to each node through majority voting after harmonizing the labels obtained for each layer.

We considered different single-layer community detection methods, including the Spectral Graph Wavelet Transform based multiscale method in [24], the Louvain Algorithm [23], as well as Spectral Clustering. The three methods obtained by cascading each single-layer method with label harmonization and majority voting are hereafter referred to as MultiLayer

<figure><img src="image_9.png" /><figcaption>Fig. 7. Application of the CPGM-CD method to the case of $M = 3$ observations of a $N = 45$ nodes graph with $P = 2, 3, 4$ communities ($p_s = 1 - p_d = 0.8$ and SNR = 40dB). We represent the communities recovered by the CPGM-CD method in terms of colors, where each color represents a specific community. The estimated communities perfectly match the true one.</figcaption></figure>

Algorithm 2 Majority voting algorithm

Input: $\Gamma^{(m)}$

Output: Community partition $\Gamma$, $P$

1: Majority Voting

2: **for** i=1:N **do**

3:

4: **for** p=1:P **do**

5: votes(p)=sum((logical($\gamma_i$==p)));

6: **end for**

7: $max_{majvote} \leftarrow max(votes)$

8: $index_{majvote} \leftarrow find(votes == max_{majvote})$

9: $index_{random} \leftarrow randomization$ of $index_{majvote}$

10: $\gamma_i \leftarrow index_{majvote}(index_{random})$

11: **end for**

SGWT (ML-SGWT) Multilayer Louvain Algorithm (ML-LA) and MultiLayer Spectral Clustering (ML-SC), respectively.

In addition, we consider a community detection method specifically introduced for multilayer networks in [21]. According to the method introduced by the authors we derive the signal on graph: $\mathbf{y} = \mathcal{H}(S)\mathbf{w}$, where $\mathbf{w}$ is a white excitation drawn from the normal distribution $N(0, I_N)$. With this state, the estimated covariance matrix $\hat{\mathbf{K}}_{yy}$ results in $\hat{\mathbf{K}}_{yy} = \frac{1}{M} \sum_{0}^{M} \mathbf{y}\mathbf{y}^T$. We consider a network process represented by the graph filter $H(L(m)) = (I - \beta L(m))^5$, with $\beta = 1/(4+4\gamma)\log N$, $\gamma = 0.3$.

We evaluate the community detection performance by computing the Normalized Mutual Information (NMI) between the detected and ground truth community labels [53]. The NMI measures the similarity between the estimated communities and the actual ones.

In the first case, we consider the observation model in (18), where the $M$ observed adjacency matrices differ due to Gaussian estimation error. We examine different Signal-to-Noise-Ratio (SNR) levels, specifically 50, 40, 30, 20, and 10 dB. To enhance statistical validity, the NMI is averaged over 500 runs. Table II presents the NMI value for the CPGM-CD method and its competitors. For the multiscale ML-SGWT method, we report results for various values of the scale parameter $s$. The CPGM-CD consistently achieves the highest NMI values compared to all the alternatives, not only at
<table><caption>TABLE II
NORMALIZED MUTUAL INFORMATION (NMI) BETWEEN THE GROUND TRUTH AND ESTIMATED COMMUNITY LABELS FOR THE CPGM-CD METHOD AND FOR STATE-OF-THE-ART COMMUNITY DETECTION METHODS, FOR A PPM GRAPH WITH *M* OBSERVATIONS OF THE GROUND-TRUTH MATRIX IN PRESENCE OF ADDITIVE NORMALLY DISTRIBUTED ESTIMATION ERROR AS IN (18) (*N*=30, *P*=2, *M*=3, *p*<sub>s</sub>=1-*p*<sub>d</sub>=0.8); THE NMI IS AVERAGED OVER 500 RUNS.</caption><thead><tr><th>NMI vs SNR</th><th>50dB</th><th>40dB</th><th>30dB</th><th>20dB</th><th>10dB</th></tr></thead><tbody><tr><td>ML-LA [23]</td><td><b>1</b></td><td>0.49</td><td>0.97</td><td>0.18</td><td>0.04</td></tr><tr><td>ML-SC [43]</td><td>0.87</td><td>0.24</td><td>0.47</td><td>0.33</td><td><b>0.07</b></td></tr><tr><td>ML-SGWT (s=1) [24]</td><td>0.03</td><td>0.02</td><td>0.02</td><td>0.03</td><td>0.03</td></tr><tr><td>ML-SGWT (s=2) [24]</td><td>0.02</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>ML-SGWT (s=3) [24]</td><td>0.02</td><td>0.05</td><td>0.09</td><td>0.03</td><td>0.04</td></tr><tr><td>ML-SGWT (s=4) [24]</td><td><b>1</b></td><td>0.02</td><td>0.13</td><td>0.06</td><td>0.04</td></tr><tr><td>ML-SGWT (s=5) [24]</td><td><b>1</b></td><td>0.86</td><td>0.62</td><td>0.17</td><td>0.05</td></tr><tr><td>ML-SGWT (s=6) [24]</td><td><b>1</b></td><td><b>1</b></td><td>0.56</td><td>0.19</td><td>0.05</td></tr><tr><td>Exact Blind [21]</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>Exact Blind (M → ∞) [21]</td><td><b>1</b></td><td>0.52</td><td>0.36</td><td>0.27</td><td>0.05</td></tr><tr><td>CPGM-CD</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.39</b></td><td><b>0.07</b></td></tr></tbody></table>

<table><caption>TABLE III
NORMALIZED MUTUAL INFORMATION (NMI) BETWEEN THE GROUND TRUTH AND ESTIMATED COMMUNITY LABELS FOR THE CPGM-CD METHOD AND FOR STATE-OF-THE-ART COMMUNITY DETECTION METHODS, FOR A PPM GRAPH AS IN (19) (*N*=30, *P*=2, *M*=3, *p*<sub>s</sub>=1-*p*<sub>d</sub>=0.8); THE NMI IS AVERAGED OVER 500 RUNS.</caption><thead><tr><th>NMI vs SNR</th><th>Noise-free</th><th>40dB</th><th>15dB</th><th>5dB</th><th>-5dB</th></tr></thead><tbody><tr><td>ML-LA [23]</td><td>0.99</td><td>0.99</td><td>0.98</td><td>0.86</td><td>0.15</td></tr><tr><td>ML-SC [43]</td><td>0.97</td><td>0.97</td><td>0.89</td><td>0.72</td><td>0.19</td></tr><tr><td>ML-SGWT (s=1) [24]</td><td>0.03</td><td>0.02</td><td>0.02</td><td>0.03</td><td>0.03</td></tr><tr><td>ML-SGWT (s=2) [24]</td><td>0.02</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>ML-SGWT (s=3) [24]</td><td>0.05</td><td>0.05</td><td>0.06</td><td>0.07</td><td>0.05</td></tr><tr><td>ML-SGWT (s=4) [24]</td><td>0.43</td><td>0.47</td><td>0.25</td><td>0.16</td><td>0.07</td></tr><tr><td>ML-SGWT (s=5) [24]</td><td>0.95</td><td>0.94</td><td>0.83</td><td>0.45</td><td>0.13</td></tr><tr><td>ML-SGWT (s=6) [24]</td><td>0.98</td><td>0.97</td><td>0.88</td><td>0.59</td><td>0.13</td></tr><tr><td>Exact Blind [21]</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>Exact Blind (M → ∞) [21]</td><td>0.94</td><td>0.94</td><td>0.93</td><td>0.79</td><td>0.15</td></tr><tr><td>CPGM-CD</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.99</b></td><td><b>0.5</b></td></tr></tbody></table>

high SNR but also at low SNR. Let us remark that for the Spectral Graph Wavelet Transform method [24], we use the un-normalized graph Laplacian due to its greater robustness to noise. The normalized graph Laplacian produces better NMI values for 50, 40, 30 dB but fails for lower SNR values. Thus, in the presence of noisy observations, the algebraic structure of the CPGM-CD effectively mitigates the error in the estimation of the adjacency matrices.

We now evaluate the performance of the CPGM-CD under the following observation model:

$$A^{(m)} = A_{GT}^{(m)} + W^{(m)}, m = 0, \dots, M-1 \quad (19)$$

where *M* ground truth adjacency matrices $A_{GT}^{(m)}$ are independently and randomly drawn under the same community partition $\Gamma$ and observed in presence of i.i.d additive Gaussian noise $W^{(m)}$, $m = 0, \dots, M-1$. Thereby, the observed adjacency matrices $A^{(m)}$ share the ground-truth community structure $\Gamma$ but they differ even at high SNR, due to sample noise. This scenario is realistic because it captures the statistical variability of real networked data, which are often collected in different experimental trials and may exhibit different network structures while sharing the same community structure. It is also challenging, as the multilayer graph constructed from the observed adjacency matrices deviates from the Cartesian Product model underlying CPGM-CD, even in the noise free case.

We report the NMI achieved by the different methods in Table III, for the case of $N = 30$, $P = 2$, $M = 3$ with 500 repetitions in order to improve the statistical validity of the simulations, at $SNR = Noise-free, 40, 15, 5, -5$. We recognize that also in this realistic and challenging scenario CPGM-CD achieves the maximum NMI on a wide range of SNR levels.

The results in Tables II,III clearly show that the CPGM-CD accurately detects the graph communities and outperforms state-of-the-art alternatives, leveraging a Cartesian Product model of the observations obtained by different estimators and/or in different trials. Let us remark that NMI values are

generally higher in Table III compared to Table II because the scenario modeled through (19) has more assorted observations and it is favorable from the point of view of community detection.

With the aim to analyze the variability of the obtained performance in function of the different parameters that can be selected to generate synthetic graphs, we now consider different scenarios. We firstly consider PPM graphs as in (18) with different number of nodes $N = 100, 150, 200$ with $P = 10, M = 5$, $p_s = 1 - p_d = 0.8$ and $SNR = 30dB$; the NMI between ground truth and estimated community labels for the different community detection methods is computed and then averaged over 100 runs. The results are shown in Fig. 8 where the case $N = 100$ is represented with a yellow line, $N = 150$ in green and $N = 200$ in light blue. From the figure, we can remark that, keeping the other parameters fixed, it is easier to correctly detect communities in a bigger network. In addition, our CPGM-CD method outperforms the state-of-the-art alternatives in the different scenarios.

Then, we consider the same graph generation scenario with PPM graphs as in (18) with different probabilities of having links between the same community $p_s = 1 - p_d = 0.7, 0.8, 0.9$ with $N = 150$ with $P = 10, M = 5$, and $SNR = 30dB$; the NMI between ground truth and estimated community labels for the different community detection methods is computed and then averaged over 100 runs. The results are presented in Fig. 9 where the case $p_s = 0.7$ is shown in yellow, $p_s = 0.8$ in green and $p_s = 0.9$ in light blue. Analyzing the results, we derive that, keeping the other parameters fixed, correctly detecting communities is easier as the intra-cluster link probability increases. Furthermore, our CPGM-CD method outperforms the state-of-the-art competitors in detecting communities in the different cases.

With the final aim to understand how the number of graph observations impacts the community detection algorithm, we consider the generation scenario with PPM graphs as in (19) varying the number of graph observations $M = 2, 5, 7$ with $N = 150$, $p_s = 1 - p_d = 0.8$ with $P = 10$ and $SNR = 1dB$. In this scenario the NMI between ground truth
<figure><img src="image_10.png" /><figcaption>Fig. 8. Normalized Mutual Information (NMI) between the ground truth and estimated community labels for the CPGM-CD method and for state-of-the-art community detection methods, for a PPM graph as in (18) with different number of nodes $N = 100, 150, 200$ with $P=10, M=5, p_s = 1 - p_d = 0.8$ and $SNR = 30dB$; the NMI is averaged over 100 runs.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Normalized Mutual Information (NMI) between the ground truth and estimated community labels for the CPGM-CD method and for state-of-the-art community detection methods, for a PPM graph as in (18) with different value of $p_s = 0.7, 0.8, 0.9$ and $N = 150, P = 10, M = 5$, and $SNR = 30dB$; the NMI is averaged over 100 runs.</figcaption></figure>

and estimated community labels for the different community detection methods is computed and then averaged over 100 runs. The results are presented in Fig. 10 where the case $M = 2$ is shown in yellow, $M = 5$ in green and $M = 7$ in light blue. We can see that, fixing the other parameters, having more graph observations generally helps the algorithms to detect communities. We also remark that the proposed CPGM-CD algorithm outperforms the state-of-the-art competitors in detecting the right communities in the different scenarios.

In addition to the theoretical foundations previously discussed, which justify the choice of the percentage of eigenvectors for the clustering stage, we now compare the NMI results using various percentages of eigenvectors. The outcomes are presented in Table IV, where four cases are reported with

<figure><img src="image_12.png" /><figcaption>Fig. 10. Normalized Mutual Information (NMI) between the ground truth and estimated community labels for the CPGM-CD method and for state-of-the-art community detection methods, for a PPM graph as in (19) with different value of layers $M = 2, 5, 7$ and $N = 150, P = 10, p_s = 1 - p_d = 0.8$ and $SNR = 1dB$; the NMI is averaged over 100 runs.</figcaption></figure>

<table><caption>TABLE IV</caption><tr><td></td><td>10%</td><td>20%</td><td>30%</td><td>50%</td></tr><tr><td>NMI</td><td>0.58</td><td>0.85</td><td>0.82</td><td>0.78</td></tr></table>

NORMALIZED MUTUAL INFORMATION (NMI) BETWEEN THE GROUND TRUTH AND ESTIMATED COMMUNITY LABELS FOR THE CPGM-CD METHOD FOR A PPM GRAPH AS IN (18) VARYING THE PERCENTAGE OF EIGENVECTORS USED FOR THE CLUSTERING (10%, 20%, 30%, 50%) WITH $N = 100, P = 10, M = 5, p_s = 1 - p_d = 0.8$ AND $SNR = 30dB$; THE NMI IS AVERAGED OVER 100 RUNS.

10%, 20%, 30%, 50% of eigenvectors associated to smallest eigenvalues. We can note that the CPGM approach remains robust regardless of the percentage of selected eigenvectors.

A few remarks are in order. In a nutshell, the performance of CPGM depends on how much the fluctuations of the observations $A^{(m)}$ impact the low-frequency subspace of the multilayer eigenvector matrix $U$. This depends on network-specific factors, namely i) the random variability of the $N \times N$ graph edges, which introduces population noise depending on the network size $N$, and ii) the ground-truth network connectivity, which affects the robustness of the multilayer low-frequency subspace. Further work is required to relate the accuracy of CPGM-CD to the ground-truth network structure (e.g. $N, P, \Gamma$ as well as to the observation model as in (18) and (19). We summarize the main trends observed in the experiments, leaving the analytical performance evaluation for future study. Firstly, for given PPM parameters $p_s, p_d$, large network structures (large $N$) result in a higher number of observed links, leading to better accuracy. Secondly, the number of estimated eigenvectors $\tilde{U}_{CD}$ used in the clustering stage depends on the product $M \cdot N_{CD}$. Therefore, a larger number of observations $M$ can be traded with a smaller number of multilayer eigenvectors $N_{CD}$. This is advantageous, as the subspace of low-frequency eigenvectors in $U$ is robust to noise [54].

Finally, when the number of communities $P$ is small relative to $N$, the community partition is more effectively revealed by
the given number of observed links.

## VI. APPLICATION ON REAL BCI EXPERIMENTS

Here, we present experimental results of our proposed community detection algorithm applied to real data recorded during motor-imagery BCI experiments. Specifically, we use open EEG data related to the BCI competition IV dataset 1 [55], involving a group of healthy subjects performing a motor imagery task without feedback. For our study, we consider motor imagery of the right hand and resting state with a total of 100 trials per condition.

This study serves as a proof-of-concept for the applicability of the proposed community detection to brain data. For simplicity, we consider data from just one subject.

We apply the CPGM-CD method in this context to identify communities associated with motor-imagery tasks.

### A. Building the multilayer network with real EEG data

The first step consists in building the multilayer network from the original EEG data. For each mental state, $N = 59$ EEG time series at 100 Hz sampling rate are obtained at the electrodes locations [55]. We compute different functional connectivity (FC) estimators on the EEG signals at nodes' pairs [9]. Since each FC estimator measures the statistical interaction between brain areas in a different way, the resulting FC graph representations differ.

In this work, we apply three classical FC estimators, i.e., the spectral coherence [56], the imaginary coherence [57] and the phase difference [58], to characterize and detect brain states during motor imagery tasks [58].

Let us consider the EEG time series collected at the graph vertices $i$ and $j$ and let us denote by $\hat{P}_i(\omega_k)$, $\hat{P}_j(\omega_k)$ and $\hat{P}_{ij}(\omega_k)$ their auto-spectra and cross-spectrum at the frequencies $\omega_k$, $k = 0, \dots, K-1$, respectively. For each nodes pair $i, j$, the spectral coherence $C_{ij}$ [56], the imaginary coherence $IC_{ij}$ [57] and the phase difference $\Phi_{ij}$ [58] are computed as the following spectral averages:

$$C_{ij} = \frac{1}{K} \sum_{k=0}^{K-1} |\hat{P}_{ij}(\omega_k)| \left( \hat{P}_i(\omega_k) \cdot \hat{P}_j(\omega_k) \right)^{-1}, \quad (20)$$

$$IC_{ij} = \frac{1}{K} \sum_{k=0}^{K-1} \left| \mathrm{Im}\{\hat{P}_{ij}(\omega_k)\} \right| \left( \hat{P}_i(\omega_k) \cdot \hat{P}_j(\omega_k) \right)^{-1}, \quad (21)$$

$$\Phi_{ij} = \frac{1}{K} \sum_{k=0}^{K-1} \left| \frac{\hat{P}_i(\omega_k)}{\hat{P}_j(\omega_k)} - \frac{\hat{P}_j(\omega_k)}{\hat{P}_i(\omega_k)} \right|. \quad (22)$$

We consider $K = 14$ frequencies $\omega_k$, corresponding to equispaced points in the $\alpha$ frequency band (8–13 Hz) involved in motor imagery tasks [59], [60].

In the top row of Fig.11 we represent the three matrices for an example of the first trial of resting state. It is possible to remark that the meaning and the scales change according to the FC estimator used for the analysis. In order to apply the Cartesian Product Graph Model, we threshold the matrices to keep a fixed percentage of the strongest connections, and rescale the values as follows. Let us remark that other methods

<figure><img src="image_13.png" /><figcaption>Fig. 11. Top line: Three FC connectivity estimators used for the model, i.e., coherence, imaginary coherence, phase difference. Bottom line: the adjacency matrix of the multilayer network, after rescaling and thresholding of each FC estimator to keep 5% of its links. In this case, we consider one trial of resting state. In the principal diagonal, we have $A^{(0)}, A^{(1)}, A^{(2)}$ derived from the FC computation with the three estimators.</figcaption></figure>

have been recently proposed to threshold the adjacency matrix [61], but we decided for the easiest and most intuitive solution.

The three observations $A^{(m)}, m=0, 1, 2$ are obtained from the different FC estimators after recasting the values in the range $[0, 1]$, with 1 corresponding to maximally correlated signals. This step is necessary in order to have coherent connectivity estimates, which corresponds to having multiple graph observations according to our Cartesian Product Graph Model. Specifically, for any $i, j, i \neq j$, we set

$$A_{ij}^{(0)} = C_{ij}, \quad A^{(1)} = 1 - IC_{ij}/M_{IC}, \quad A_{ij}^{(2)} = 1 - \Phi_{ij}/M_{\Phi},$$

being $M_{IC} = \max_{ij} IC_{ij}$ and $M_{\Phi} = \max_{ij} \Phi_{ij}$; besides, $A_{ii}^{(m)} = 0, m=0, \dots, M-1, i=0, \dots, N-1$.

Fig.11 shows in the bottom line, the multilayer adjacency matrix $A$ after thresholding and rescaling the original FC estimators, and we can see the three observed matrices $A^{(m)}, m = 0, 1, 2$ on the principal block diagonal and the upper and lower block diagonals are the identity matrices, corresponding to the interlayer connections.

### B. Community detection of the multilayer network with real EEG data

In this subsection, we apply our proposed community detection algorithm to real EEG data previously described. In order to investigate the community organization during different human mental states, we consider EEG data of one subject performing two different tasks, that are: motor-imagery of the right hand and resting state. We remark that this study represents a proof of concept to show the potential of our method to identify the community structure on brain data in absence of ground truth community labels. Interestingly, the analysis can be extended to a large number of subjects to evaluate the impact of the estimated community structure in a
classification system such as a Brain Computer interface; this is left for further study.

In order to understand the functioning of the method, we present in Fig. 12 the result of the community detection algorithm on resting state. In particular, we consider FC estimators computed in $\alpha$ frequency band and we replicate the analysis for the 100 available trials of resting state. In order to have a compact but representative visualization, we derive the most representative trial. To this aim, we compute the NMI for each pair of trials, since it measures the similarity between communities obtained at different observations. Then, for each trial, we sum all the derived NMIs and we take the one associated to the largest sum. This operation identifies the most representative trial because it is the most similar to all the others. Results are reported in terms of colors, which means that each node of the network is coloured according to the community identified by the algorithm. In this case, $P = 5$. In Fig. 12 the link width represents the stability between the three networks: indeed, the link is widest when it is present (after thresholding) on the three networks and it is absent if there is no link in any graph. For what concerns the connections, we recognize that the stable links between nodes, represented by thick gray lines, are captured by the multilayer mode, and the CPGM-CD detects communities, represented by nodes of the same colors, that reflect such stable links. More specifically, we recognize links between frontal nodes (in the blue community), links between parietal nodes (green community) and links in left, right and central areas, respectively in orange, red and light blue.

<figure><img src="image_14.png" /><figcaption>Fig. 12. Result of the community detection algorithm in terms of communities represented on graph for resting state in $\alpha$ band for all the 100 trials. The colour of the node represents the community to whom it belongs. We represent only the 2% of strongest links at each layer. The thickness of the links represent the stability across the different layers, i.e. the thicker the link, the more stable it is across the three estimators considered.</figcaption></figure>

In order to investigate the organization of the communities on the brain scalp, we represent them in a different representation in the top line of Fig. 13. In the three figures of panel

(a), we have different views of the brain scalp where each node is coloured according to the associated communities. The identified clusters appear symmetrically distributed in the brain and they geometrically follow the brain organization: we have a community in the frontal area, represented in blue, a community localized in the central electrodes, in light blue, then we have two symmetric communities in the two hemispheres, i.e., coloured in orange and in red, and a community that contains nodes in parietal and temporal electrodes (represented in yellow).

We now compare the cluster found during resting state with another condition, i.e., motor imagery of the right hand. Results of the analysis related to the most representative trial are reported in the panel (b) of Fig. 13. The output of the algorithm differs from the resting state condition. Indeed, here we remark the presence of the community in the left hemisphere that widens including all the left hemisphere in the contralateral motor area, that is supposed to be involved in the motor imagery task, until the center and the beginning of the right hemisphere. Then we have a community in the frontal electrodes, two communities in the parietal lobe and one community, represented in orange, located in the right sensori-motor cortex. Our findings show that during motor imagery i) the symmetry of the brain communities is disrupted in the sensori-motor area, and ii) the community associated to the contralateral motor area widens.

### C. Statistical analysis of community detection-based features with real EEG data

In this subsection, we apply our CPGM-CD method to the previously described EEG data to characterize mental states associated with motor imagery in different subjects. This analysis proves that the organization of brain communities serves as a distinguishing feature between these states. To construct a matrix that captures the community structure while remaining independent of trial-dependent labeling, we represent the $p$-th community by an $N$-dimensional indicator vector, $b_p$, where the $i$-th element is equal to 1 if the $i$-th node belongs to the $p$-th community, and 0 otherwise. We then compute the matrix $\Psi$, defined as follows:

$$ \Psi = \sum_{p=0}^{P-1} b_p \cdot b_p^T \qquad (23) $$

With these positions, $\Psi$ is an $N \times N$ matrix where each binary element $\Psi_{ij}$ indicates whether the two nodes $i$ and $j$ belong to the same community (1) or not (0).

The advantage of computing the matrix $\Psi$ is that it remains independent of the community labeling order, which may vary across different experiments.

We quantify the impact of community-based features for distinguishing between two mental states of a BCI user during a motor imagery task, by: i) collecting the matrix $\Psi$ over different trials, and ii) performing statistical tests to compare the distribution of its elements, $\Psi_{ij}$, under two conditions-motor imagery of the left hand and resting. The discriminative ability of the community-based matrix $\Psi$ is assessed by the
<figure><img src="image_15.png" /><figcaption>Fig. 13. Output of community detection algorithm on real EEG data, where the colors of each node correspond to the resulting community. In the top line, we have results of EEG data recorded during motor imagery of the left hand, and in the bottom line we have those related to motor imagery of the right hand. For this analysis, we previously selected *P* = 4 communities and we bandpassed signals in the alpha frequency band.</figcaption></figure>

$t$-values, which indicate whether the element $\Psi_{ij}$ shows a statistical distance between the two conditions. In this context, the $t$-value quantifies the difference in the $\Psi$ matrix, measured across 100 repeated trials, for a fixed frequency band and a specific subject.

In Fig. 14(a)-(c), we present the results of our analysis in terms of $t$-values with a significance value of $\rho = 0.05$ for three subjects, within $\alpha$ frequency band. The highlighted pairs of $i, j$-th nodes are characterized by differences in their corresponding $\Psi_{ij}$ values- indicating whether the nodes belong to the same community or not- across the two conditions.

From the results, we observe that all subjects exhibit a significant difference in motor imagery of the left hand and the resting state. This suggests that certain nodes belong to the same community in one mental state, while in the other state, the community organization shifts, causing previously linked nodes to no longer be connected. Notably, the areas of significant changes are the sensori-motor cortex, and the motor area contralateral to the imagined movement, specifically the right hemisphere. This indicates that connectivity changes and differences in community organizations occur in the right hemisphere, which is expected to be involved in this particular motor imagery task.

Our study serves as a proof of concept for the applicability and potential of a community detection approach on brain networks. It is important to note that fine-tuning the frequency

band and electrode selection, along with appropriate pre-processing steps, is essential to evaluate the practical applicability of the proposed method.

Overall, our approach effectively captures the variability in human brain reconfiguration during motor imagery experiments. Therefore, brain clusters could be used as an additional feature for controlling a real Brain-Computer Interface, providing novel insights beyond those offered by traditional power spectral density or functional connectivity measures.

## VII. CONCLUSIONS

This paper introduces an original model for representing multiple observations using on a multilayer network that can be described as the Cartesian Product Graph between the ground truth graph- the actual graph underlying the possibly noisy observations- and a path graph, which has as many nodes as the number of layers. The key contribution of this work is the introduction of a novel approach for identifying community structures in networked data by estimating optimized eigenvectors of the multilayer graph for community detection purpose. Our method has been shown to outperform state-of-the-art competitors in detecting community structures on synthetic data. Furthermore, we apply the proposed framework to real EEG data collected during motor imagery- based BCI experiments, as a proof of concept with the aim to show its potential to reveal functionally relevant communities in the
<figure><img src="image_16.png" /><figcaption>Fig. 14. Results of statistical tests between motor imagery of the left hand vs. rest across 100 trials from three different subjects (a),(b),(c).</figcaption></figure>

human brain. Further research and parameter tuning is needed to extend the analysis to a bigger dataset and other tasks to properly assess the applicability of the method on real EEG data during real-time BCIs.

Overall, the proposed method shows promise for modeling the brain EEG network and detecting community structures during cognitive tasks. It is also applicable in cases where different estimates of a graph adjacency matrix are available, paving the way for further generalization of the observation model.

# APPENDIX A
## COMPUTATIONAL COMPLEXITY

The computational complexity of our method depends on the estimation of the first $N_{CD}$ eigenvectors. It can be approximated to $O(N_{Edges} \cdot N_{CD})$, where $N_{Edges}$ is the total number of edges and $N_{CD}$ is the number of eigenvectors that have to be estimated. In our case, the total number of edges $N_{Edges}$ depends both on the total number of intralayer links $N_e$ and the interlayer links $N \cdot (M-1)$.

Specifically, the total number of edges is intralayer links $N_e \cdot M$, because each layer has $N_e$ edges and, according to the proposed Cartesian Product Graph Model, the total number of intelayer connections is $N \cdot (M-1)$, where $N$ is the number of nodes of each graph and $M$ is the number of layers. The computational cost $C$ corresponds to $C = N_e \cdot N \cdot M(M-1)$. Let us underline here that the graph is sparse, $N_e \ll N \cdot M$.

<figure><img src="image_17.png" /><figcaption>Fig. 15. Example of harmonization result before majority voting. In this case, the parameters are: $N = 25$, $P = 2$, $M = 2$, $N_{runs} = 1$, $prob = 0.9$. In the top line, we report the communities obtained with Louvain procedure with different colors before the harmonization. In the bottom line, we have the communities (identified with colors) after the harmonization step.</figcaption></figure>

# APPENDIX B
## LABEL HARMONIZATION

The centroids of the clusters found at the first iteration are computed and the subsequent iterations the new centroids are compared with the stored ones and the name of the labels are possibly changed according to the smallest geometrical distance between novel and stored labels. In Fig. 15, we graphically represent the effect of the label harmonization. Indeed, in the top line, we report the clusters identified with the Louvain method, with $N = 25$, $P = 2$, $M = 2$, $N_{runs} = 1$, $prob = 0.9$ at the different layers. It is possible to remark that the clusters are exactly the same but the labels are shifted, since the community represented in blue at the first layer correspond to the red at the second layer. In the bottom line, we represent the cluster after the harmonization stage.

# ACKNOWLEDGMENT

FD acknowledges support from the European Research Council (ERC), the European Union's through Horizon 2020 research and innovation programme under Grant 864729

# REFERENCES

[1] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, "Graph signal processing: Overview, challenges, and applications," *Proceedings of the IEEE*, vol. 106, no. 5, pp. 808–828, 2018.

[2] A. Sandryhaila and J. M. Moura, "Discrete signal processing on graphs," *IEEE transactions on signal processing*, vol. 61, no. 7, pp. 1644–1656, 2013.

[3] S. Fortunato, "Community detection in graphs," *Physics reports*, vol. 486, no. 3-5, pp. 75–174, 2010.

[4] A. Lancichinetti and S. Fortunato, "Community detection algorithms: a comparative analysis," *Physical review E*, vol. 80, no. 5, p. 056117, 2009.

[5] W. Huang, T. A. Bolton, J. D. Medaglia, D. S. Bassett, A. Ribeiro, and D. Van De Ville, "A graph signal processing perspective on functional brain imaging," *Proceedings of the IEEE*, vol. 106, no. 5, pp. 868–885, 2018.
[6] C.-M. Ting, S. B. Samdin, M. Tang, and H. Ombao, "Detecting dynamic community structure in functional brain networks across individuals: A multilayer approach," *IEEE Transactions on Medical Imaging*, vol. 40, no. 2, pp. 468-480, 2020.

[7] M. Ménoret, N. Farrugia, B. Pasdeloup, and V. Gripon, "Evaluating graph signal processing for neuroimaging through classification and dimensionality reduction," in *2017 IEEE Global Conference on Signal and Information Processing (GlobalSIP)*, pp. 618-622, IEEE, 2017.

[8] S. Itani and D. Thanou, "A graph signal processing framework for the classification of temporal brain data," in *2020 28th European Signal Processing Conference (EUSIPCO)*, pp. 1180-1184, IEEE, 2021.

[9] A. M. Bastos and J.-M. Schoffelen, "A tutorial review of functional connectivity analysis methods and their interpretational pitfalls," *Frontiers in systems neuroscience*, vol. 9, p. 175, 2016.

[10] J. Guillon, Y. Attal, O. Colliot, V. La Corte, B. Dubois, D. Schwartz, M. Chavez, and F. de Vico Fallani, "Loss of brain inter-frequency hubs in alzheimer's disease," *Scientific reports*, vol. 7, no. 1, pp. 1-13, 2017.

[11] C. Presigny and F. D. V. Fallani, "Colloquium: Multiscale modeling of brain network organization," *Reviews of Modern Physics*, vol. 94, no. 3, p. 031002, 2022.

[12] S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, "The structure and dynamics of multilayer networks," *Physics reports*, vol. 544, no. 1, pp. 1-122, 2014.

[13] T. Cattai, G. Scarano, M.-C. Corsi, F. D. V. Fallani, and S. Colonnese, "Eeg as signal on graph: a multilayer network model for bci applications," in *2022 30th European Signal Processing Conference (EUSIPCO)*, pp. 932-936, IEEE, 2022.

[14] G. Ortiz-Jiménez, M. Coutino, S. P. Chepuri, and G. Leus, "Sampling and reconstruction of signals on product graphs," in *2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP)*, pp. 713-717, IEEE, 2018.

[15] H. Sayama, "Estimation of laplacian spectra of direct and strong product graphs," *Discrete Applied Mathematics*, vol. 205, pp. 160-170, 2016.

[16] J. Kim and J.-G. Lee, "Community detection in multi-layer graphs: A survey," *ACM SIGMOD Record*, vol. 44, no. 3, pp. 37-48, 2015.

[17] 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, 2014.

[18] W. Tang, Z. Lu, and I. S. Dhillon, "Clustering with multiple graphs," in *2009 Ninth IEEE International Conference on Data Mining*, pp. 1016-1021, IEEE, 2009.

[19] X. Dong, P. Frossard, P. Vandergheynst, and N. Nefedov, "Clustering with multi-layer graphs: A spectral perspective," *IEEE Transactions on Signal Processing*, vol. 60, no. 11, pp. 5820-5831, 2012.

[20] Z. Zeng, J. Wang, L. Zhou, and G. Karypis, "Coherent closed quasi-clique discovery from large dense graph databases," in *Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining*, pp. 797-802, 2006.

[21] T. M. Roddenberry, M. T. Schaub, H.-T. Wai, and S. Segarra, "Exact blind community detection from signals on multiple graphs," *IEEE Transactions on Signal Processing*, vol. 68, pp. 5016-5030, 2020.

[22] P. De Meo, E. Ferrara, G. Fiumara, and A. Provetti, "Generalized louvain method for community detection in large networks," in *2011 11th international conference on intelligent systems design and applications*, pp. 88-93, IEEE, 2011.

[23] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, "Fast unfolding of communities in large networks," *Journal of statistical mechanics: theory and experiment*, vol. 2008, no. 10, p. P10008, 2008.

[24] N. Tremblay and P. Borgnat, "Graph wavelets for multiscale community mining," *IEEE Transactions on Signal Processing*, vol. 62, no. 20, pp. 5227-5239, 2014.

[25] M. T. Schaub, S. Segarra, and J. N. Tsitsiklis, "Blind identification of stochastic block models from dynamical observations," *SIAM Journal on Mathematics of Data Science*, vol. 2, no. 2, pp. 335-367, 2020.

[26] L. Peel, J.-C. Delvenne, and R. Lambiotte, "Multiscale mixing patterns in networks," *Proceedings of the National Academy of Sciences*, vol. 115, no. 16, pp. 4057-4062, 2018.

[27] R. F. Betzel, M. A. Bertolero, E. M. Gordon, C. Gratton, N. U. Dosenbach, and D. S. Bassett, "The community structure of functional brain networks exhibits scale-specific patterns of inter- and intra-subject variability," *NeuroImage*, vol. 202, p. 115990, 2019.

[28] M. E. Newman, "Assortative mixing in networks," *Physical review letters*, vol. 89, no. 20, p. 208701, 2002.

[29] B. Karrer and M. E. Newman, "Stochastic blockmodels and community structure in networks," *Physical review E*, vol. 83, no. 1, p. 016107, 2011.

[30] P. W. Holland, K. B. Laskey, and S. Leinhardt, "Stochastic blockmodels: First steps," *Social networks*, vol. 5, no. 2, pp. 109-137, 1983.

[31] T. Cattai, G. Scarano, M.-C. Corsi, D. S. Bassett, F. D. V. Fallani, and S. Colonnese, "Improving j-divergence of brain connectivity states by graph laplacian denoising," *IEEE transactions on Signal and Information Processing over Networks*, vol. 7, pp. 493-508, 2021.

[32] A. Milanese, J. Sun, and T. Nishikawa, "Approximating spectral impact of structural perturbations in large networks," *Physical Review E*, vol. 81, no. 4, p. 046112, 2010.

[33] A. Michaël, *Eigenvalues and eigenvectors of large matrices under random perturbations*. PhD thesis, Université Sorbonne Paris Cité, 2018.

[34] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains," *IEEE signal processing magazine*, vol. 30, no. 3, pp. 83-98, 2013.

[35] B. Ricaud, P. Borgnat, N. Tremblay, P. Gonçalves, and P. Vandergheynst, "Fourier could be a data scientist: From graph fourier transform to signal processing on graphs," *Comptes Rendus Physique*, vol. 20, no. 5, pp. 474-488, 2019. Fourier and the science of today / Fourier et la science d'aujourd'hui.

[36] F. Xia, K. Sun, S. Yu, A. Aziz, L. Wan, S. Pan, and H. Liu, "Graph learning: A survey," *IEEE Transactions on Artificial Intelligence*, vol. 2, no. 2, pp. 109-127, 2021.

[37] L. Qiao, L. Zhang, S. Chen, and D. Shen, "Data-driven graph construction and graph learning: A review," *Neurocomputing*, vol. 312, pp. 336-351, 2018.

[38] G. Ortiz-Jiménez, M. Coutino, S. P. Chepuri, and G. Leus, "Sampling and reconstruction of signals on product graphs," in *2018 IEEE Global Conference on Signal and Information Processing*, pp. 713-717, 2018.

[39] S. White and P. Smyth, "A spectral clustering approach to finding communities in graphs," in *Proceedings of the 2005 SIAM international conference on data mining*, pp. 274-285, SIAM, 2005.

[40] L. Su, W. Wang, and Y. Zhang, "Strong consistency of spectral clustering for stochastic block models," *IEEE Transactions on Information Theory*, vol. 66, no. 1, pp. 324-338, 2019.

[41] U. Von Luxburg, "A tutorial on spectral clustering," *Statistics and computing*, vol. 17, pp. 395-416, 2007.

[42] J. Shi and J. Malik, "Normalized cuts and image segmentation," *IEEE Transactions on pattern analysis and machine intelligence*, vol. 22, no. 8, pp. 888-905, 2000.

[43] A. Ng, M. Jordan, and Y. Weiss, "On spectral clustering: Analysis and an algorithm," *Advances in neural information processing systems*, vol. 14, 2001.

[44] A. Lancichinetti and S. Fortunato, "Consensus clustering in complex networks," *Scientific reports*, vol. 2, no. 1, p. 336, 2012.

[45] K. R. Shahapure and C. Nicholas, "Cluster quality analysis using silhouette score," in *2020 IEEE 7th international conference on data science and advanced analytics (DSAA)*, pp. 747-748, IEEE, 2020.

[46] T. Caliński and J. Harabasz, "A dendrite method for cluster analysis," *Communications in Statistics-theory and Methods*, vol. 3, no. 1, pp. 1-27, 1974.

[47] X. Wang and Y. Xu, "An improved index for clustering validation based on silhouette index and calinski-harabasz index," in *IOP Conference Series: Materials Science and Engineering*, vol. 569, p. 052024, IOP Publishing, 2019.

[48] T. Cattai, C. Caporali, M.-C. Corsi, and S. Colonnese, "Introducing the modularity graph: an application to brain functional networks," in *2024 32nd European Signal Processing Conference (EUSIPCO)*, pp. 1611-1615, IEEE, 2024.

[49] H.-W. Shen and X.-Q. Cheng, "Spectral methods for the detection of network community structure: a comparative analysis," *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2010, no. 10, p. P10020, 2010.

[50] F. de Vico Fallani, J. Richiardi, M. Chavez, and S. Achard, "Graph analysis of functional brain networks: practical issues in translational neuroscience," *Philosophical Transactions of the Royal Society B: Biological Sciences*, vol. 369, no. 1653, p. 20130521, 2014.

[51] E. Kolaczyk, "Analysis of network data: methods and models," 2009.

[52] M. Browne and A. Shapiro, "The asymptotic covariance matrix of sample correlation coefficients under general conditions," *Linear Algebra and its Applications*, vol. 82, pp. 169-176, 1986.
[53] N. X. Vinh, J. Epps, and J. Bailey, "Information theoretic measures for clusterings comparison: is a correction for chance necessary?," in *Proceedings of the 26th annual international conference on machine learning*, pp. 1073–1080, 2009.

[54] T. Cattai, G. Scarano, M.-C. Corsi, D. S. Bassett, F. De Vico Fallani, and S. Colonnese, "Improving j-divergence of brain connectivity states by graph laplacian denoising," *IEEE Transactions on Signal and Information Processing over Networks*, vol. 7, pp. 493–508, 2021.

[55] M. Tangermann, K.-R. Müller, A. Aertsen, N. Birbaumer, C. Braun, C. Brunner, R. Leeb, C. Mehring, K. J. Miller, G. Mueller-Putz, et al., "Review of the bci competition iv," *Frontiers in neuroscience*, vol. 6, p. 55, 2012.

[56] G. C. Carter, "Coherence and time delay estimation," *Proceedings of the IEEE*, vol. 75, no. 2, pp. 236–255, 1987.

[57] G. Nolte, O. Bai, L. Wheaton, Z. Mari, S. Vorbach, and M. Hallett, "Identifying true brain interaction from eeg data using the imaginary part of coherency," *Clinical neurophysiology*, vol. 115, no. 10, pp. 2292–2307, 2004.

[58] T. Cattai, S. Colonnese, M.-C. Corsi, D. S. Bassett, G. Scarano, and F. D. V. Fallani, "Phase/amplitude synchronization of brain signals during motor imagery bci tasks," *IEEE Transactions on Neural Systems and Rehabilitation Engineering*, vol. 29, pp. 1168–1177, 2021.

[59] G. Pfurtscheller and C. Neuper, "Motor imagery and direct brain-computer communication," *Proceedings of the IEEE*, vol. 89, no. 7, pp. 1123–1134, 2001.

[60] G. Pfurtscheller, C. Brunner, A. Schlögl, and F. L. Da Silva, "Mu rhythm (de) synchronization and eeg single-trial classification of different motor imagery tasks," *NeuroImage*, vol. 31, no. 1, pp. 153–159, 2006.

[61] F. De Vico Fallani, V. Latora, and M. Chavez, "A topological criterion for filtering information in complex brain networks," *PLoS computational biology*, vol. 13, no. 1, p. e1005305, 2017.