# Recovering Static and Time-Varying Communities Using Persistent Edges

Konstantin Avrachenkov, Maximilien Dreveton, and Lasse Leskelä

**Abstract**—This article focuses on spectral methods for recovering communities in temporal networks. In the case of fixed communities, spectral clustering on the simple time-aggregated graph (i.e., the weighted graph formed by the sum of the interactions over all temporal snapshots) does not always produce satisfying results. To utilise information carried by temporal correlations, we propose to employ different weights on freshly appearing and persistent edges. We show that spectral clustering on such weighted graphs can be explained as a relaxation of the maximum likelihood estimator of an extension of the degree-corrected stochastic block model with Markov interactions. We also study the setting of evolving communities, for which we use the prediction at time $t-1$ as an oracle for inferring the community labels at time $t$. We demonstrate the accuracy of the proposed methods on synthetic and real data sets.

**Index Terms**—Graph clustering, temporal networks, spectral methods, stochastic block model.

## I. INTRODUCTION

COMPLEX networks are commonly used to describe and analyze interactions between entities. A natural problem arising in complex network analysis consists of identifying meaningful structures within the network. Community recovery, i.e., partitioning the set of nodes of a network into communities based on some common properties of the vertices, is now a well-established area [1]. In many situations, interactions between node pairs vary over time, and classical graph-based models are replaced by temporal network models [2]. Temporal networks naturally arise in various situations, such as communication or face-to-face interactions between individuals. In brain networks, interactions between neurons are also time-varying. Finally, ecological networks such as food webs experience seasonal changes. The longitudinal dimension of data raises new challenges to traditional clustering algorithms. Moreover, the community structure of temporal networks can either be static

or time-varying, and different clustering algorithms need to be used in each case. For example, community structures in a brain network are likely to be static, while the communities in a social network can vary.

If communities are static, then each additional snapshot makes clustering easier. Furthermore, when snapshots are temporally uncorrelated, simple time-aggregation of data does not lose any information and recovers communities optimally [3]. Nonetheless, for time-correlated models, this simple temporal aggregation of the data might lose important features such as temporal patterns. As such, [4] describes information-theoretic recovery criteria for Markov edge evolution dynamics, while [5], [6] study a spectral algorithm based on the squared adjacency matrix.

Previous research on community detection in temporal networks has also focused on evolving communities, for example by generalizing belief-propagation methods [7], developing variational EM algorithms [8] and introducing new spectral methods [9], [10], [11] or modularity-based methods [12], [13]. We refer to [14], [15] for an overview. All of the aforementioned works focus on evolving communities for which the interactions between nodes are re-sampled at every time step. One can then treat each layer independently by applying static community detection and smoothing the community predictions. A more challenging situation is when the edges at a given time step also depend on the edges at the previous time step. It is argued in [16] that link persistence makes the inference of communities harder.

In this work, we first study the case of static community memberships. We introduce a temporal extension of the degree-corrected stochastic block model [17], [18], in which the community structure is fixed and the interactions between node pairs follow a Markov evolution which only depends on the community structure and the degree correction parameters. We show that maximum likelihood estimation reduces to the maximisation of a regularized modularity of a time-aggregated graph, in the limit of a large number of snapshots and sparse interactions. This weighted graph is not given by a simple sum of the adjacency matrices over all snapshots. Instead, we show that the maximum likelihood estimator distinguishes, at each snapshot, the newly formed edges from the persistent ones. Indeed, the information contained in these two types of temporal edges is different. Thus, we propose a time-aggregation corresponding to the sum of all snapshots of the signal coming from the newly formed and persistent edges, but with different weights. This leads to a weighted graph that better takes into account the temporal nature of the signal than the vanilla time-aggregation.

Manuscript received 24 April 2023; revised 30 October 2023; accepted 22 November 2023. Date of publication 30 November 2023; date of current version 23 February 2024. This work was supported in part by the Project of Inria - Nokia Bell Labs "Distributed Learning and Control for Network Analysis" and in part by COSTNET EU Cost Action CA15109. Recommended for acceptance by Dr. Luoyi Fu. (Corresponding author: Konstantin Avrachenkov.)

Konstantin Avrachenkov is with Inria Sophia Antipolis, 06902 Valbonne, France (e-mail: k.avrachenkov@inria.fr).

Maximilien Dreveton is with the School of Computer and Communication Sciences, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland (e-mail: maximilien.dreveton@epfl.ch).

Lasse Leskelä is with the Department of Mathematics and Systems Analysis, Aalto University, 02150 Espoo, Finland (e-mail: lasse.leskela@aalto.fi).

Digital Object Identifier 10.1109/TNSE.2023.3337281
A continuous relaxation of this maximum likelihood estimator then leads to a normalized spectral clustering algorithm. We finally extend this study to higher-order temporal correlations.

We also study evolving communities with temporal correlation in the edge formation. In this case, we propose to recover the community at a given time step by using the prediction of the previous time step as a noisy oracle. This naturally leads to a noisy semi-supervised inference problem [15], which we tackle using spectral methods.

The article is structured as follows. Section II studies the static communities, while Section III focuses on evolving communities. Numerical experiments are presented in Section IV and the conclusions and future research directions are given in Section V.

This article is an extension of [19]. Compared to the conference article, the section on evolving communities is entirely new. Moreover, the spectral algorithm for fixed communities is modified (in order not to require any hyper-parameters besides the number of communities). The extension to higher-order Markov chains is also new. Finally, we present experiments on a wider collection of data sets.

## II. STATIC COMMUNITIES

This section studies the recovery of static communities in a temporal network. We introduce in Section II-A a temporal stochastic block model where dynamic interaction patterns between node pairs follow a Markov chain. We render this model more versatile by adding degree correction parameters, describing the tendency of each node to start new interactions. We show that in some cases the likelihood of this model is approximated by a regularized modularity of a time-aggregated graph. This time-aggregated graph involves a trade-off between new edges and persistent edges. A continuous relaxation reduces the regularized modularity maximisation to a normalized spectral clustering.

### A. Degree-Corrected Temporal Network Model

A degree-corrected temporal stochastic block model with N nodes, K blocks and T snapshots is a probability distribution

$$ \mathbb{P}(A | \sigma, F, \theta) = \prod_{1 \le i < j \le N} F_{\sigma_i \sigma_j}^{\theta_i \theta_j} (A_{ij}^1, \dots, A_{ij}^T) \quad (1) $$

of a symmetric adjacency tensor $A \in \{0, 1\}^{N \times N \times T}$ with zero diagonal entries, defined by $A_{ij}^t = 1$ if nodes i and j interact at time t, and $A_{ij}^t = 0$ otherwise. Moreover, $\sigma = (\sigma_1, \dots, \sigma_N)$ is a community assignment with $\sigma_i \in \{1, \dots, K\}$ indicating the community of node i, $F = (F_{kl}^{\eta \theta})$ is a collection of probability distributions on $\{0, 1\}^T$ indexed by community labels k, l $\in [K] = \{1, \dots, K\}$ and numbers $\eta, \theta \ge 0$ (representing the degree correction parameters), and $\theta = (\theta_1, \dots, \theta_N)$ is a vector of node-specific degree correction parameters $\theta_i \ge 0$.

In the following, we will restrict ourselves to homogeneous models with Markov edge dynamics, in which the nodes' static

community memberships are sampled independently and uniformly at random from $[K]$, and for $x \in \{0, 1\}^T$ we have

$$ F_{\sigma_i \sigma_j}^{\theta_i \theta_j}(x) = \begin{cases} \mu_{x_1}^{\theta_i \theta_j} \prod_{t=2}^{T} P_{x_{t-1}, x_t} & \text{if } \sigma_i = \sigma_j, \\ \nu_{x_1}^{\theta_i \theta_j} \prod_{t=2}^{T} Q_{x_{t-1}, x_t}^{\theta_i \theta_j} & \text{otherwise,} \end{cases} \quad (2) $$

with initial distributions

$$ \mu^{\theta_i \theta_j} = \begin{pmatrix} 1 - \theta_i \theta_j \mu_1 \\ \theta_i \theta_j \mu_1 \end{pmatrix}, \\ \nu^{\theta_i \theta_j} = \begin{pmatrix} 1 - \theta_i \theta_j \nu_1 \\ \theta_i \theta_j \nu_1 \end{pmatrix}, \quad (3) $$

and transition probability matrices

$$ P^{\theta_i \theta_j} = \begin{pmatrix} 1 - \theta_i \theta_j P_{01} & \theta_i \theta_j P_{01} \\ 1 - P_{11} & P_{11} \end{pmatrix}, \\ Q^{\theta_i \theta_j} = \begin{pmatrix} 1 - \theta_i \theta_j Q_{01} & \theta_i \theta_j Q_{01} \\ 1 - Q_{11} & Q_{11} \end{pmatrix}. \quad (4) $$

The parameters $\theta_i$ account for the fact that some nodes might be more inclined than others to start new connections, similarly to the degree-corrected block model of [18]. More precisely, $\theta_i \theta_j \mu_1$ (resp., $\theta_i \theta_j \nu_1$) is the probability of an interaction between nodes i and j belonging to the same (resp., different) cluster at the first snapshot. The first snapshot is therefore an instance of the degree-corrected block model of [18]. Moreover, the interaction between i and j is a Markov chain whose transition matrix is $P^{\theta_i \theta_j}$ (resp., $Q^{\theta_i \theta_j}$) if the two nodes belong to the same (resp., different) cluster. To obtain tractable analytical expressions for the model likelihood, we do not add degree correction parameters in front of $P_{11}$ and $Q_{11}$. This corresponds to assuming that popularity does not affect the duration of interactions. Moreover, we assume that $\max_{i,j}\{\theta_i \theta_j \delta\} \le 1$, where $\delta = \max\{\mu_1, \nu_1, P_{01}, Q_{01}\}$. Finally, we normalise the degree correction parameters so that $\frac{1}{|C_k|} \sum_{i \in C_k} \theta_i = 1$ for all i, where $C_k = \{i : \sigma(i) = k\}$ indicates community k.

### B. Maximum Likelihood Estimator

The following proposition presents a concise formulation of the model likelihood, with proof presented in Appendix A.

**Proposition 2.1:** A maximum likelihood estimator for the Markov block model defined by (1)–(2) is any community assignment $\sigma \in [K]^N$ that maximises

$$ \sum_{\substack{i,j \\ \sigma_i = \sigma_j}} \left\{ A_{ij}^1 (\rho_1^{\theta_i \theta_j} - \rho_0^{\theta_i \theta_j}) + \rho_0^{\theta_i \theta_j} + (A_{ij}^1 - A_{ij}^T) \ell_{10}^{\theta_i \theta_j} \right\} \\ + \sum_{\substack{i,j \\ \sigma_i = \sigma_j}} \sum_{t=2}^{T} \left\{ (\ell_{01}^{\theta_i \theta_j} + \ell_{10}^{\theta_i \theta_j}) (A_{ij}^t - A_{ij}^{t-1} A_{ij}^t) \right. \\ \left. + \ell_{11}^{\theta_i \theta_j} A_{ij}^{t-1} A_{ij}^t - \log \frac{Q_{00}^{\theta_i \theta_j}}{P_{00}^{\theta_i \theta_j}} \right\} $$
where $\rho_a^{\theta_i \theta_j} = \log \frac{\mu_a^{\theta_i \theta_j}}{\nu_a^{\theta_i \theta_j}}$ and $\ell_{ab}^{\theta_i \theta_j} = \log \frac{P_{ab}^{\theta_i \theta_j}}{Q_{ab}^{\theta_i \theta_j}} - \log \frac{P_{00}^{\theta_i \theta_j}}{Q_{00}^{\theta_i \theta_j}}$.

The MLE derived in Proposition 2.1 is more complex than summing all snapshots independently. In particular, the terms $A_{ij}^{t-1} A_{ij}^t$ account for persistent edges over two consecutive snapshots. Denote by $A_{\text{pers}}^t = A^{t-1} \odot A^t$ the entrywise product of adjacency matrices $A^{t-1}$ and $A^t$. Then $A_{\text{pers}}^t$ is the adjacency matrix of the graph containing the persistent edges between $t-1$ and $t$, and $A_{\text{new}}^t = A^t - A_{\text{pers}}^t$ corresponds to the graph containing the edges freshly appearing at time $t$.

Assuming that the number of snapshots $T$ is large, we can ignore the boundary terms, and the MLE expressed in Proposition 2.1 reduces to maximising

$$ \sum_{t=2}^{T} \sum_{\substack{i,j: \\ \sigma_i=\sigma_j}} \left\{ \left( \ell_{01}^{\theta_i \theta_j} + \ell_{10}^{\theta_i \theta_j} \right) (A_{ij}^t - A_{ij}^{t-1} A_{ij}^t) \right. \\ \left. + \ell_{11}^{\theta_i \theta_j} A_{ij}^{t-1} A_{ij}^t - \log \frac{Q_{00}^{\theta_i \theta_j}}{P_{00}^{\theta_i \theta_j}} \right\}. $$

By utilising (3)–(4), we can further simplify it to express this as a modularity. Recall that given a weighted graph $W$, a partition $\sigma$, and a resolution parameter $\gamma$, the regularized modularity is defined as [20], [21]

$$ \mathcal{M}(W, \sigma, \gamma) = \sum_{i,j} \delta(\sigma_i, \sigma_j) \left( W_{ij} - \gamma \frac{d_i d_j}{2m} \right), $$

where $d_i = \sum_j W_{ij}$ and $m = \sum_i d_i$. Hence, suppose that $P^{\theta_i \theta_j}$ and $Q^{\theta_i \theta_j}$ are nondegenerate, and $\mu^{\theta_i \theta_j}$ (resp. $\nu^{\theta_i \theta_j}$) is the stationary distribution of $P^{\theta_i \theta_j}$ (resp. $Q^{\theta_i \theta_j}$). In a sparse setting, $P_{01}$ and $Q_{01}$ are small, and after a Taylor expansion (see Section A-B for detailed derivations) the previous expression is approximately equal to $\mathcal{M}(W, \sigma, \gamma)$, where $W$ is defined by

$$ W = \sum_{t=2}^{T} \left( \alpha A_{\text{new}}^{t} + \beta A_{\text{pers}}^{t} \right), \quad (5) $$

with

$$ \begin{align} \alpha &= \log \frac{P_{01}}{Q_{01}} + \log \frac{1 - P_{11}}{1 - Q_{11}}, \nonumber \\ \beta &= \log \frac{P_{11}}{Q_{11}}, \tag{6} \end{align} $$

$$ \gamma = (P_{01} - Q_{01}) \frac{\alpha(\mu_1 + (K-1)\nu_1) + (\beta - \alpha)(\mu_1 P_{11} + (K-1)\nu_1 Q_{11})}{K}. $$

### C. Comparison With Previous Work on Static Communities

The correspondence between maximum likelihood estimation and modularity maximisation is known in static block models [22]. Analogously to the single-layer case, the modularity of a temporal network, with possibly time-dependent community structure, was previously defined by Mucha et al. [12] and Pamfil et al. [13] as

$$ \sum_{t=1}^{T} \mathcal{M}(A^t, \sigma^t, \gamma_t) + \sum_{t=1}^{T} \sum_{s \neq t} \sum_{i} \omega_i^{st} \delta(\sigma_i^s, \sigma_i^t) \quad (7) $$

where $\gamma_t$ is the resolution parameter for layer $t$, $\sigma_i^t$ is the community membership of node $i$ at time step $t$, and $w_i^{st}$ denotes a coupling between time instants $s$ and $t$. For a static community structure, the second term in (7) is irrelevant. When the resolution is constant over time, the relevant term in (7) can be written as

$$ \sum_{t=1}^{T} \mathcal{M}(A^t, \sigma, \gamma) = \mathcal{M}(A^{\text{agg}}, \sigma, \gamma), $$

where $A^{\text{agg}} = \sum_{t=1}^{T} A^t$ is the weighted adjacency matrix of the time-aggregated data. In contrast, the matrix $W$ in (5) involves a trade-off between new edges and persistent edges. We notice that $W \propto A^{\text{agg}}$ only if $\alpha = \beta$. This is the case when freshly appearing and persistent edges are equally important, and only in this case using the simple time-aggregated adjacency matrix provides optimal inference.

### D. Temporal Spectral Clustering Combining New and Persistent Edges

Following our analysis in Section II-B, the community prediction is given by

$$ \hat{\sigma} = \arg \max_{\sigma \in [K]^N} \mathcal{M}(W, \sigma, \gamma), $$

where $W$ is defined in (5) and $\gamma$ is a proper resolution parameter. This optimisation problem is NP-complete in general [23] but can be approximately solved by continuous relaxation. We can choose the relaxation so that the optimisation problem reduces to normalized spectral clustering on the weighted graph $W$ (we refer to [24] and to Section A-C for detailed computations). We note that in order to compute the normalized Laplacian of $W$, we should assume $\alpha, \beta \ge 0$, which we observe in all our experiments.

Since $\alpha$ and $\beta$ are likely to be unknown in practical situations, we propose to estimate them as follows. For $a, b \in \{0, 1\}$, the empirical probability of observing a transition $a \to b$ in the interaction pattern between nodes $i$ and $j$ is

$$ \frac{n_{ab}(i, j)}{n_a(i, j)}, $$

where $n_{ab}(i, j) = \sum_{t=2}^{T} 1(A_{ij}^{t-1} = a, A_{ij}^t = b)$ is the observed number of transitions $a \to b$ between $i$ and $j$, and $n_a(i, j) = \sum_b n_{ab}(i, j)$. Thus, if $\hat{\sigma}$ is an estimator of the community labelling, then

$$ \begin{align} \hat{P}_{ab} &= \frac{1}{|\{i, j : \hat{\sigma}_i = \hat{\sigma}_j\}|} \sum_{i,j: \hat{\sigma}_i = \hat{\sigma}_j} \frac{n_{ab}(i, j)}{n_a(i, j)}, \nonumber \\ \hat{Q}_{ab} &= \frac{1}{|\{i, j : \hat{\sigma}_i \neq \hat{\sigma}_j\}|} \sum_{i,j: \hat{\sigma}_i \neq \hat{\sigma}_j} \frac{n_{ab}(i, j)}{n_a(i, j)}, \tag{8} \end{align} $$

are estimators of $P$ and $Q$ from which one can compute estimates of $\alpha$ and $\beta$. These observations motivate the following algorithm. We start with $\alpha = \beta = 1$, and we let $\hat{\sigma}$ be the output of spectral clustering on $W = \sum_{t=2}^{T} \alpha A_{\text{new}}^t + \beta A_{\text{pers}}^t$. Using $\hat{\sigma}$ we compute estimates $\hat{P}$ and $\hat{Q}$, which lead to a new value of $\alpha$ and $\beta$. We
**Algorithm 1:** Spectral Clustering for Temporal Networks
With Markov Edge Dynamics and Static Node Labelling.

**Input:** Adjacency matrices $A^1, \dots, A^T$, number of clusters $K$, stopping parameter $\epsilon$

**Output:** Predicted membership matrix $\hat{\sigma} \in [K]^N$

**Process:**

Let $\alpha = \beta = 1$ and stopping = False;

while stopping = False do

• Let $W = \sum_{t=2}^{T} \alpha A_{\text{new}}^{t} + \beta A_{\text{pers}}^{t}$ where $A_{\text{new}}^{t} = A^{t} - A^{t-1} \odot A^{t}$ and $A_{\text{pers}}^{t} = A^{t-1} \odot A^{t}$;

• $\hat{\sigma} \leftarrow \text{SPECTRAL CLUSTERING}(W, K)$;

• Let $\hat{P}$ and $\hat{Q}$ as defined by (2.8);

• Let $\tilde{\alpha} = \log \frac{\hat{P}_{01}}{\hat{Q}_{01}} + \log \frac{\hat{P}_{10}}{\hat{Q}_{10}}$ and $\tilde{\beta} = \log \frac{\hat{P}_{11}}{\hat{Q}_{11}}$.

if $\|(\alpha, \beta) - (\tilde{\alpha}, \tilde{\beta})\|_2 < \epsilon$ then

else

$$ \lfloor \alpha = \tilde{\alpha} \text{ and } \beta = \tilde{\beta}. $$

**Return:** $\hat{\sigma}$

iterate the procedure until some convergence criterion is met. This leads to Algorithm 1.

### E. Extension to Higher-Order Markov Chains

Interactions with a longer memory can be modelled using higher-order Markov chains [25], [26]. We demonstrate below how Algorithm 1 can be extended to a long-memory context. Suppose the interactions are Markov chains of order $M$. Namely, the intra- and inter-block interaction connectivity functions, defined previously by (2), are now given by

$$ F_{\sigma_i \sigma_j}(x) = \begin{cases} \mu(x_{1:M}) \prod_{t=M+1}^{T} P(x_t | s_t) & \text{if } \sigma_i = \sigma_j, \\ \nu(x_{1:M}) \prod_{t=M+1}^{T} Q(x_t | s_t) & \text{otherwise,} \end{cases} $$

where $s_t = (x_{t-M}, \dots, x_{t-1})$. Following a similar reasoning, we propose to apply spectral clustering on the weighted matrix $W$ defined by

$$ W_{ij} = \sum_{t=M+1}^{T} \sum_{\substack{s \in \{0,1\}^M \\ b \in \{0,1\}}} \alpha_{sb} 1(A_{ij}^{(t-M):(t-1)} = s, A_{ij}^{t} = b), $$

with $\alpha_{sb} = \log \frac{P(b|s)}{Q(b|s)}$ and $A_{ij}^{(t-M):(t-1)} = (A_{ij}^{t-M}, \dots, A_{ij}^{t-1})$. Moreover, an estimator of $P$ is given by

$$ \hat{P}(b|s) = \frac{1}{|\{i,j : \hat{\sigma}_i = \hat{\sigma}_j\}|} \sum_{\hat{\sigma}_i = \hat{\sigma}_j} \frac{n_{sb}(i,j)}{n_s(i,j)}, $$

$$ \hat{Q}(b|s) = \frac{1}{|\{i,j : \hat{\sigma}_i \neq \hat{\sigma}_j\}|} \sum_{\hat{\sigma}_i \neq \hat{\sigma}_j} \frac{n_{sb}(i,j)}{n_s(i,j)}, \quad (9) $$

where $n_{sb}(i,j) = \sum_{t=M}^{T} 1(A_{ij}^{(t-M):(t-1)} = s, A_{ij}^{t} = b)$ and $n_s(i,j) = \sum_{b \in \{0,1\}} n_{sb}(i,j)$. This can be formalised as Algorithm 2.

**Algorithm 2:** Spectral Clustering for Temporal Networks
With Higher-Order Markov Edge Dynamics.

**Input:** Adjacency matrices $A^1, \dots, A^T$, number of clusters $K$, Markov order $M$, stopping parameter $\epsilon$

**Output:** Predicted membership matrix $\hat{\sigma} \in [K]^N$

**Process:**

Let $\alpha = (\alpha_{sb})_{s \in \{0,1\}^M, b \in \{0,1\}}$ and initialise $\alpha_{sb} = 1$ for all $s, b$. Let stopping = False;

while stopping = False do

• $W_{ij} := \sum_{t=2}^{T} \sum_{s,b} \alpha_{sb} 1(A_{ij}^{(t-M):(t-1)} = s, A_{ij}^{t} = b)$;

• $\hat{\sigma} \leftarrow \text{SPECTRAL CLUSTERING}(W, K)$;

• Let $\hat{P}$ and $\hat{Q}$ as defined by (2.9);

• Let $\tilde{\alpha}_{sb} := \log \frac{\hat{P}}{\hat{Q}}(b|s)$ for all $s \in \{0,1\}^M$ and $b \in \{0,1\}$.

if $\|\alpha - \tilde{\alpha}\|_2 < \epsilon$ then

else
$$ \lfloor \alpha = \tilde{\alpha}. $$

**Return:** $\hat{\sigma}$

### III. EVOLVING COMMUNITIES

In this section, we consider a population of $N$ nodes partitioned into $K$ time-evolving communities. At time $t$, we denote by $\sigma_i^t \in [K]$ the community membership of node $i$ and by $A_{ij}^t \in \{0,1\}$ the observed interaction between nodes $i$ and $j$. We investigate methods of recovering the community structure, denoted by $(\sigma^1, \dots, \sigma^T)$ where $\sigma^t \in [K]^N$, from an observed adjacency tensor $A = (A_{ij}^t)$.

#### A. Model Description

Similarly to several articles on dynamic SBM [7], [8], [16], we first assume that each node community labels $\sigma_i^{1:T} \in [K]^T$ is a Markov chain of length $T$ with initial probability $\alpha$ and transition probability matrix $\pi$. Hence,

$$ \mathbb{P}(\sigma^{1:T}) = \prod_{i=1}^{N} \alpha(\sigma_i^1) \prod_{t=2}^{T} \pi(\sigma_i^{t-1}, \sigma_i^t). \quad (10) $$

For simplicity, we will assume that the initial labels and the transitions are uniform, that is

$$ \alpha = \frac{1}{K} 1_K \quad \text{and} \quad \pi = \eta I_K + \frac{1-\eta}{K} 1_K 1_K^T. $$

In other words, a node keeps its label with probability $\eta \in [0, 1]$ and chooses a label uniformly at random with probability $1 - \eta$.

We then assume that the interaction between two nodes $i$ and $j$ is a Markov chain depending only on the community labelling.<fnref n="1" /> Hence,

$$ \mathbb{P}(A|\sigma) = \prod_{1 \le i < j \le N} \mathbb{P}(A_{ij}^1 | \sigma_i^1, \sigma_j^1) \prod_{t=2}^{T} \mathbb{P}(A_{ij}^t | A_{ij}^{t-1}, \sigma_i^t, \sigma_j^t). $$

<footnote n="1">¹For simplicity, in this section we do not add any degree correction parameters.</footnote>
We consider a homogeneous model in which the initial distribution is given by

$$ \mathbb{P}(A_{ij}^1 | \sigma_i^1, \sigma_j^1) = \begin{cases} \mu(A_{ij}^1), & \text{if } \sigma_i^1 = \sigma_j^1, \\ \nu(A_{ij}^1), & \text{otherwise,} \end{cases} $$

and the transition probabilities are

$$ \mathbb{P}(A_{ij}^t = b | A_{ij}^{t-1} = a, \sigma_i^t, \sigma_j^t) = \begin{cases} P_{ab} & \text{if } \sigma_i^t = \sigma_j^t, \\ Q_{ab} & \text{otherwise.} \end{cases} $$

## B. Community Inference

The inference of the temporal community structure is complicated because interactions observed in a distant past are less informative for inferring the current community membership of a node. To mitigate the bias caused by outdated information, we propose a two-stage online algorithm:

1) at time $t = 1$, we use a static community detection algorithm to output $\hat{\sigma}^1$, a prediction of the initial node labels $\sigma^1$ from the observation of the first snapshot $A^1$;

2) at time $t > 1$, we will use the observation of the first $t$ snapshots $A^1, \dots, A^t$ as well as the previous predictions $\hat{\sigma}^1, \dots, \hat{\sigma}^{t-1}$. This will be treated as a semi-supervised learning problem, where the prediction $\hat{\sigma}^{t-1}$ is seen as a noisy oracle for the true node labelling $\sigma^t$.

From the Markov structure, the prediction at time $t > 1$ reduces to predicting $\sigma^t$ using only the network at time $t-1$ and $t$ and the previous prediction $\hat{\sigma}^{t-1}$. This can be interpreted as a noisy semi-supervised problem, where the previous prediction $\hat{\sigma}^{t-1}$ plays the role of an oracle for the node labels at time $t$. This oracle is biased, as it bears two kinds of potential mistakes. Firstly, $\hat{\sigma}^{t-1}$ is not necessarily exactly equal to the perfect community labelling $\sigma^{t-1}$. Furthermore, since the node labels vary through time, $\sigma^{t-1}$ does not precisely correspond to $\sigma^t$. Assume that the network data $A$ and community labels $\sigma$ come from the model described in Section III-A. We denote the oracle bias at time $t$ by

$$ \rho_t = \mathbb{P}(\sigma_i^t \neq \hat{\sigma}_i^{t-1}), \quad (11) $$

and we will make a simplifying assumption that it does not depend on node $i$.

The MAP estimator for the online learning problem is

$$ \begin{aligned} \hat{\sigma}^t &= \underset{z \in [K]^N}{\arg \max} \mathbb{P} (z | A^t, A^{t-1}, \hat{\sigma}^{t-1}) \\ &= \underset{z \in [K]^N}{\arg \max} \mathbb{P} (A^t | A^{t-1}, z, \hat{\sigma}^{t-1}) \mathbb{P} (z | A^{t-1}, \hat{\sigma}^{t-1}) \end{aligned} $$

by Bayes' rule.

Since $\mathbb{P}(A^t | A^{t-1}, z, \hat{\sigma}^{t-1}) = \mathbb{P}(A^t | A^{t-1}, z)$, then by proceeding similarly to the proof of Proposition 2.1, the log-likelihood term $\log \mathbb{P}(A^t | A^{t-1}, z)$ can be rewritten as

$$ \frac{1}{2} \sum_{\substack{i,j \\ z_i=z_j}} \left\{ \begin{aligned} & \ell_{01}^{\theta_i \theta_j} (A_{ij}^t - A_{ij}^{t-1} A_{ij}^{t-1}) + \ell_{10}^{\theta_i \theta_j} (A_{ij}^{t-1} - A_{ij}^{t-1} A_{ij}^t) \\ & + \ell_{11}^{\theta_i \theta_j} A_{ij}^{t-1} A_{ij}^t - \log \frac{Q_{00}^{\theta_i \theta_j}}{P_{00}^{\theta_i \theta_j}} \end{aligned} \right\}. $$

The oracle information is equal to

$$ \begin{aligned} \mathbb{P}(z | \hat{\sigma}^{t-1}) &= \prod_{i=1}^{N} \frac{\mathbb{P}(\hat{\sigma}_i^{t-1} | z_i)}{\mathbb{P}(\hat{\sigma}_i^{t-1})} \mathbb{P}(z_i) \\ &= \left( \frac{\rho_t}{1 - \rho_t} \right)^{|\{i \in [N]: z_i \neq \hat{\sigma}_i^{t-1}\}|} (1 - \rho_t)^N K^{-N}, \end{aligned} $$

where we used the uniformity of the node labels and where $\rho_t$ is defined in (11).

The MAP estimator at time $t$ is therefore any labelling $\sigma$ that maximises

$$ \sum_{\substack{i,j \\ z_i=z_j}} \left\{ \ell_{01}^{\theta_i \theta_j} (A_{ij}^t - A_{ij}^{t-1} A_{ij}^{t-1}) + \ell_{10}^{\theta_i \theta_j} (A_{ij}^{t-1} - A_{ij}^{t-1} A_{ij}^t) \right. \\ \left. \qquad + \ell_{11}^{\theta_i \theta_j} A_{ij}^{t-1} A_{ij}^t - \log \frac{Q_{00}^{\theta_i \theta_j}}{P_{00}^{\theta_i \theta_j}} \right\} + 2\lambda_t \sum_{i=1}^{N} 1(z_i = s_i), $$

where $\ell_{ab}^{\theta_i \theta_j} = \log \frac{P_{ab}^{\theta_i \theta_j}}{P_{ab}^{\theta_i \theta_j}} - \log \frac{P_{00}^{\theta_i \theta_j}}{P_{00}^{\theta_i \theta_j}}$ and $\lambda_t = \log \frac{1-\rho_t}{\rho_t}$, and where $s_i = \sigma_i^{t-1}$.

## C. Continuous Relaxation

Denote by $A_{\text{pers}}^t = A^{t-1} \odot A^t$ the persistent edges, by $A_{\text{new}} = A^t - A_{\text{pers}}^t$ the freshly formed edges, and by $A_{\text{old}} = A^{t-1} - A_{\text{pers}}^t$ the disappearing edges between time $t-1$ and $t$. Denote by $\Sigma \in \{0, 1\}^{N \times K}$ (resp., by $S$) the one-hot encoding of $\sigma \in [K]^N$ (resp., of $s$). In other words, we have $\Sigma_{ik} = 1(\sigma_i = k)$ and $S_{ij} = 1(s_i = k)$. We observe that $\sum_i 1(\sigma_i = s_i) = \Sigma^T S$. Thus, using a Taylor expansion as in Section II-A, we can approximate the MAP estimator at time $t$ as the maximisation over $\Sigma \in \{0, 1\}^{N \times K}$ of

$$ \mathrm{Tr} (\Sigma^T (W^t - \tau 1_N 1_N^T) \Sigma + 2\lambda \Sigma^T S), \quad (12) $$

where $W^t = \alpha_{01} A_{\text{new}}^t + \alpha_{10} A_{\text{old}}^t + \alpha_{11} A_{\text{pers}}^t$ with $\alpha_{ab} = \log \frac{P_{ab}^{\theta_i \theta_j}}{P_{ab}^{\theta_i \theta_j}}$ and $\tau = P_{01} - Q_{01}$ is a resolution parameter, and $\Sigma$ verifies the constraints $\sum_{k \in [K]} \Sigma_{ik} = 1$ for all $i \in [N]$.

We perform a continuous relaxation (mirroring what is commonly done for spectral methods [24]) and solve instead

$$ \begin{array}{ll} \underset{X \in \mathbb{R}^{N \times K}}{\arg \max} & \mathrm{Tr} (X^T W_{\tau} X + 2\lambda X^T S), \quad (13) \\ & \sum_{i=1}^{N} \kappa_i X_{ik}^2 = \sum_{i=1}^{N} \kappa_i \\ & \text{where } \kappa = (\kappa_1, \dots, \kappa_N) \text{ is a vector of positive entries. In what follows, we choose to constrain } X \text{ to verify } (X^T DX)_{kk} = 1 \text{ by taking } \kappa_i = \frac{d_i}{2m}.^2 \text{ This leads to the optimisation problem} \end{array} $$

$$ \underset{X \in \mathbb{R}^{N \times K}}{\arg \max} \quad \mathrm{Tr} (X^T W_{\tau} X + 2\lambda X^T S). \quad (14) \\ \forall k \in [K] : (X^T DX)_{kk} = 1 $$

<footnote n="2">²Another possible relaxation is to constrain $X$ such that $X^T DX = I_K$, so that $D^{-1/2}X$ belongs to a Stiefel manifold. The minimisation problem obtained with this (stronger) constraint is similar to an unbalanced orthogonal Procrustes problem, for which necessary or sufficient conditions for the local and/or global minimiser exist, but no closed-form expression of the solution is available [27], [28].</footnote>
**Algorithm 3:** Online Clustering of Time-Varying Communities.

**Input:** Observed graph sequence

$A^{1:T} = (A^1, \dots, A^T)$; number of communities $K$; static graph clustering algorithm algo; parameters $\alpha_{01}, \alpha_{10}, \alpha_{11}$ and $\lambda_1, \dots, \lambda_T$

**Output:** Node labelling $\hat{\sigma}^{1:T}$

**Initialise:** Compute $\hat{\sigma}^1 \leftarrow \text{algo}(A^1)$.

1 **for** $t = 2, \dots, T$ **do**

2 Compute $W = \alpha_{01}A_{\text{new}}^t + \alpha_{10}A_{\text{old}}^t + \alpha_{11}A_{\text{pers.}}^t$

3 For all $k \in [K]$, let $\gamma_k^*$ be the smallest solution of Equation (3.8).

4 Compute $\hat{X}_{k}$ as the solution of Equation (3.7).

5 Let $\hat{\sigma}_i^t = \text{K-means}(\hat{X}_{i\cdot})$ for all $i \in [N]$.

In order to solve (14), let us introduce the eigenvalue decomposition of $-D^{-1/2}W_{\tau}D^{-1/2}$, written as

$$-D^{-1/2}W_{\tau}D^{-1/2} = Q\Delta Q^T, \quad (15)$$

where $\Delta = \text{diag}(\delta_1, \dots, \delta_N)$ with $\delta_1 \le \dots \le \delta_N$ and $Q^TQ = I_N$.

The following Proposition, whose proof is deferred to Appendix B, gives an expression for the solution of (14).

**Proposition 3.1:** Let $B = \lambda Q^T D^{-1/2} S$ where $Q$ is defined in (15). The solution $\hat{X} = (\hat{X}_{1\cdot}, \dots, \hat{X}_{K\cdot})$ of the optimisation problem (14) satisfies for all $k \in [K]$

$$(-W_{\tau} - \gamma_k^* D) \hat{X}_{k\cdot} = \lambda S_{k\cdot}, \quad (16)$$

where $\gamma_k^*$ is the smallest solution of

$$\sum_{i=1}^{N} \left( \frac{B_{ik}}{\delta_i - \gamma} \right)^2 - 1 = 0. \quad (17)$$

The expression of the solution of (14) given in Proposition 3.1 leads to Algorithm 3.

## D. Comparison With Previous Work on Evolving Communities

It is often natural to assume that the pairwise interactions evolve over time, but the temporal changes in connectivity patterns occur gradually. Therefore, many works in the literature propose dynamic extensions of the stochastic block model, where both community memberships and connections can vary through time.

To model a smooth evolution of communities across time, the majority of works describe changes in the memberships via a Markov evolution. Moreover, previous works almost always suppose that the interaction at time $t$ depends only on the block structure at time $t$. For example, Yang et al. [29] assume that the community membership of each node forms a Markov chain independent of the memberships of other nodes. Xu and Hero [30] extend this model and allow for both the community memberships and connectivity parameters to vary. This renders the model unidentifiable. This issue was raised by Matias and Miele [8], who established identifiability conditions

for models with both community memberships and connectivity parameters varying over time. Pensky and Zhang [31] assume that at most $s$ nodes can switch their memberships between two consecutive time instants, and that the connection probabilities vary smoothly over time.

All of the aforementioned works suppose that the interactions at a given snapshot $t$ are conditionally independent of all past snapshots given the current community memberships. Namely, it is assumed that the temporal interactions are resampled at each time instant (i.e., the interaction $A_{ij}^t$ between nodes $i$ and $j$ at time $t$ depends only on the community labels $\sigma_i^t$ and $\sigma_j^t$ at time $t$, and not on the previous interactions $A_{ij}^{t-1}$ between $i$ and $j$). Barucca et al. [16] propose a model similar to ours but assume that the edge persistence is the same for intra- and inter-community interactions. More precisely, using the notations of Section III-A, [16] supposes that

$$\mathbb{P} (A_{ij}^t = b | A_{ij}^{t-1} = a, \sigma_i^t, \sigma_j^t) = \begin{cases} P_{ab}, & \text{if } \sigma_i^t = \sigma_j^t, \\ Q_{ab}, & \text{otherwise,} \end{cases}$$

where the transition matrices $P$ and $Q$ can be written as

$$\begin{aligned} P &= \xi_1 \mathbb{I}_2 \mathbb{I}_2 + (1 - \xi) \tilde{P}, \\ Q &= \xi_1 \mathbb{I}_2 \mathbb{I}_2 + (1 - \xi) \tilde{Q}, \end{aligned} \quad (18)$$

and where $\tilde{P}$ and $\tilde{Q}$ are given transition matrices and $\xi \in [0, 1]$ is an edge persistence parameter. Such a decomposition implicitly supposes that the Markov chain is positively autocorrelated<fnref n="3" /> ($\text{Cov}(X_0, X_1) \ge 0$) for a stationary Markov chain ($X_t$)). While positive autocorrelation can be a reasonable assumption for social networks, it might not be suitable for other situations (for example, in biological networks spiking phenomena might occur).

To generalise the model even further, in addition to time-varying interactions and community memberships, the Markov transition matrices describing the persistence of edges may also evolve over time. An instance of such a model was discussed by Xu [32], by assuming a Gaussian evolution for log-transformed transition matrices, and experimenting with an extended Kalman filter for inferring community memberships and model parameters. Further exploration of such more general Markov models appears an interesting direction for future research.

## IV. NUMERICAL EXPERIMENTS

### A. Static Communities

We compare the performances of temporal spectral clustering data sets with static communities. In this section, *improved time-aggregation* refers to Algorithm 1, while *simple time-aggregation* refers to spectral clustering on the time-aggregated graph $\sum_t A^t$.

<footnote n="3">³This comes from the following Lemma. A 2-by-2 binary stochastic matrix $P \neq I_2$ admits a representation $P = \xi I_2 + (1 - \xi) \begin{pmatrix} \pi_0 & \pi_1 \\ \pi_1 & \pi_0 \end{pmatrix}$ for some $\xi \in [0, 1]$ and some probability distribution $\pi = (\pi_0, \pi_1)$ iff both eigenvalues of $P$ are non-negative iff $\text{Cov}(X_0, X_1) \ge 0$ for a stationary Markov chain $(X_0, X_1, \dots)$ on $\{0, 1\}$ with transition matrix $P$.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1. Accuracy of Algorithm 1 with $\alpha = 1$ and different $\beta$, on a degree-corrected temporal SBM with 300 nodes and $K = 3$ blocks (with a uniform prior), and a stationary Markov edge evolution $\mu_1 = 0.02$, $\nu_1 = 0.01$, $P_{11} = 0.7$ and $Q_{11} = 0.2$, for a different generation of the degree correction parameters $\theta$. The results are averaged over 25 synthetic graphs, and error bars show the standard error.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2. Comparison of accuracy obtained on the SocioPatterns data sets by using spectral clustering with simple time-aggregation versus improved time-aggregation (Algorithm 1).</figcaption></figure>

In our experiments, we set $\epsilon = 0.1$ and a maximum number of iterations = 10. We observed empirically that Algorithm 1 frequently required less than 5 iterations to converge. The only cases when Algorithm 1 did not stop within 10 iterations were the cases when the time-aggregated graphs were too sparse (this may happen when the time horizon $T$ is not enough long).

1) Synthetic Data Sets: We plot in Fig. 1 the averaged adjusted rand-index (ARI) obtained on 25 realisations of stochastic block models with Markov edge dynamics for various degree-correction parameters $\theta$. More precisely:

* $\theta_i = 1$ for all $i \in [N]$ in Fig. 1(a);

* in Fig. 1(b), the $\theta_i$'s are generated $\theta_i$ according to $|\mathcal{N}(0, \sigma^2)| + 1 - \sigma\sqrt{2/\pi}$ where $|\mathcal{N}(0, \sigma^2)|$ denotes the

absolute value of a normal random variable with mean 0 and variance $\sigma^2$, where $\sigma = 0.25$;

* in Fig. 1(c), the $\theta_i$'s are generated from a Pareto distribution with density function $f(x) = \frac{am^a}{x^{a+1}}1(x \ge m)$ with $a = 2.5$ and $m = 2/3$.

Note that the sampling of the $\theta_i$'s enforces $E\theta_i = 1$ in all settings.

While the simple time aggregation works well, it is striking to notice that Algorithm 1 provides better results. For example, when the degree-corrected parameters are generated according to a Gaussian distribution, Fig. 1(b) shows that the simple time-aggregation requires around 100 snapshots to get a perfect ARI, while Algorithm 1 only requires 30 snapshots.
<figure><img src="image_3.png" /><figcaption>Fig. 3. Comparison of the ARI obtained on workplace 2015 and primary school data sets by using temporal spectral clustering with higher order correlations (Algorithm 2).</figcaption></figure>

<table><caption>TABLE I
DIMENSIONS OF SIX DATA SETS FROM SOCIOPATTERN LIBRARY</caption><thead><tr><th>Data set name</th><th>N</th><th>T</th><th>K</th></tr></thead><tbody><tr><td>High school 2011</td><td>118</td><td>5391</td><td>3</td></tr><tr><td>High school 2012</td><td>180</td><td>11273</td><td>5</td></tr><tr><td>High school 2013</td><td>327</td><td>7375</td><td>9</td></tr><tr><td>Workplace 2013</td><td>88</td><td>7027</td><td>4</td></tr><tr><td>Workplace 2015</td><td>161</td><td>17508</td><td>5</td></tr><tr><td>Primary school</td><td>232</td><td>3100</td><td>10</td></tr></tbody></table>

2) *Real World Data:* We investigate data sets of face-to-face contacts between individuals, for which ground-truth communities are known. The data sets are available on SocioPatterns website.<fnref n="4" /> Three data sets were collected in a high school [33], [34], one in a primary school [35] and two in a work environment [36], [37]. For each data set, nodes correspond to individuals (students or workers), time-varying network edges correspond to close proximity interactions, measured by wearable sensors able to sense proximity, and communities correspond to school classes or work departments. The data set dimensions are given in Table I.

We plot in Fig. 2 the accuracy obtained by Algorithm 1 on each of these data sets, compared to the spectral clustering done on the time aggregated graph $\sum_{t=1}^{T} A^t$. We observe a clear improvement in the accuracy of most data sets. This highlights the importance of weighting differently the persistent edges and the freshly appearing edges. Additionally, we would like to emphasise that Algorithm 1 estimates the values of $\alpha$ and $\beta$ using the network data. Only two data sets (workplace 2015 and primary school) show no improvements. Nonetheless, we demonstrate in Fig. 3 that we can improve the accuracy of these two data sets by taking into account higher-order correlations.

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

(a) $\eta = 1$

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

(b) $\eta = 0.85$

<caption>Fig. 4. ARI obtained by various clustering algorithms on time-varying Markov Block Models with 300 nodes and $K = 2$ blocks (with a uniform prior), and a stationary Markov edge evolution $\mu_1 = 0.1, \nu_1 = 0.05, P_{11} = 0.6$ and $Q_{11} = 0.2$. The results are averaged over 25 synthetic graphs, and error bars show the standard error. We ran Algorithm 3 (online-ssl) with parameters $\alpha_{01} = 1, \alpha_{10} = 0$ and $\alpha_{11} = 2$), and we compare with weighted SC (Algorithm 1 with $\alpha = 1, \beta = 2$) and an algorithm performing spectral clustering on each snapshot individually (individual SC).</caption>

<footnote n="4">⁴[Online]. Available: https://www.sociopatterns.org/</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. ARI of Algorithm 3 with $\alpha_{01} = 1$, $\alpha_{10} = 0$ and $\alpha_{11} = 2$ for various choices of $\lambda$. Simulations are performed on time-varying Markov Block Models with $n = 200$, $K = 2$, $\mu_1 = 0.2$, $\nu_1 = 0.1$, $P_{11} = 0.7$, $Q_{11} = 0.3$ and $\eta = 0.8$. The results are averaged over 25 synthetic graphs, and error bars show the standard error.</figcaption></figure>

## B. Evolving Communities

In this section, we evaluate the performance of Algorithm 3 on synthetic data sets.

We compare in Fig. 4 the averaged accuracy obtained by Algorithm 3 with Algorithm 1 (spectral clustering with persistent edges) and an algorithm performing spectral clustering on each snapshot individually. In particular, we observe that when $\eta = 1$ (i.e., static community structure), Algorithm 1 is extremely efficient, as expected. Since it takes into account all previous snapshots, it in particular outperforms Algorithm 3. On the contrary, when $\eta \neq 1$, the lagging problem arises, and Algorithm 1 ends up with a very poor accuracy after a few snapshots. On the contrary, Algorithm 3 keeps a very high accuracy over all snapshots.

In Fig. 4, we choose $\lambda_t$ to be constant with respect to time and equal to 0.5, while Fig. 5 explores other constant possible values of $\lambda$. We observe that when $\lambda$ is in the range $[0.1, 1]$, the algorithm provides good accuracy. On the other hand, when $\lambda$ becomes too large, Algorithm 3 gives too much importance to the oracle, and the accuracy becomes worse. In practice, the choice of the parameters $\lambda_t$ could be optimised from the data, e.g. based on $\eta$ or on the transition matrices $P$ and $Q$. Moreover, it would be intuitive to increase $\lambda_t$ with $t$, as the confidence in the oracle is higher when more temporal data is available. We leave these research ideas for future work.

Finally, we compare in Fig. 6 the performance of Algorithm 3 with the VEM algorithm of [8].<fnref n="5" /> The VEM algorithm of [8] is tailored for models with time-varying communities but where the interactions are resampled at every time step. Thus, we observe that its performance is excellent when the links are resampled at every time step (Fig. 6(a)), but the performance decreases when the time-persistency increases (Fig. 6(b) and (c)). In contrast,

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

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

<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison of ARIs obtained by Algorithm 3 and the VEM algorithm of [8] on time-varying Markov Block Models with 250 nodes and $K = 2$ blocks (with a uniform prior), $\eta = 0.75$, and a stationary Markov edge evolution $\mu_1 = 0.2$, $\nu_1 = 0.1$, $Q_{11} = 0.1$ and for various $P_{11}$. Fig. 6(a) is a setting where links are re-sampled at each time step, while Fig. 6(b) and (c) models medium and high intra-community link persistency. The results are averaged over 25 synthetic graphs, and error bars show the standard error.</figcaption></figure>

Algorithm 3 performs better than VEM when the correlation in the edge dynamics is strong.

## V. CONCLUSION

In this article, we derived spectral algorithms for clustering temporal networks, both in the case of static and evolving

<footnote n="5">⁵To implement this VEM algorithm, we used the R package *dynsbm* provided by the authors of [8] and at http://lbbe.univ-lyon1.fr/dynsbm.</footnote>
communities. For static communities, the Markov chain assumption for the edge formation implies a difference in the information carried by persistent and newly formed edges. This difference is lost by performing a simple time aggregation that indifferently sums all the interactions. Instead, we proposed to weigh differently the persistent and the newly formed edges. In the evolving community setting, we highlighted the possibility of using a past prediction as a semi-supervised oracle for the prediction at a current snapshot. Overall, the derived spectral algorithms for clustering temporal networks provide a valuable tool for inferring community structure in complex systems, and the proposed weighing methods for persistent and newly formed edges offer a refined approach to analyse dynamic networks.

One direction for future research is online/adaptive tuning of hyperparameters such as estimates of the bias of the noisy semi-supervised oracle.

Several approximations and relaxations were made in deriving a weighted similarity matrix used for spectral clustering and for tuning its weights. Analysing the accuracy of these approximations, and proving that the proposed algorithms terminate in a finite time remain topics for future research.

## ACKNOWLEDGMENT

We thank anonymous reviewers for helping to improve the presentation of the article.

## APPENDIX A
PROOFS FOR SECTION II

### A. MLE Computation

**Proof of Proposition 2.1:** By the temporal Markov property, the log-likelihood of the model can be written as $\log \mathbb{P}(A|\sigma, \theta) = \log \mathbb{P}(A^1|\sigma, \theta) + \sum_{t=2}^{T} \log \mathbb{P}(A^t|A^{t-1}, \sigma, \theta)$. By denoting $\rho_a^{\theta_i\theta_j} = \log \frac{\mu_a^{\theta_i\theta_j}}{\nu_a^{\theta_i\theta_j}}$, we find that

$$ 
\begin{aligned}
\log \mathbb{P}(A^1 | \sigma, \theta) \\
&= \frac{1}{2} \sum_{i,j} \sum_a \delta(A_{ij}^1, a) \left( \delta(\sigma_i, \sigma_j) \rho_a^{\theta_i\theta_j} + \log \nu_a^{\theta_i\theta_j} \right) \\
&= \frac{1}{2} \sum_{i,j} \delta(\sigma_i, \sigma_j) \sum_a \delta(A_{ij}^1, a) \rho_a^{\theta_i\theta_j} + c_1(A),
\end{aligned}
$$

where $c_1(A) = \frac{1}{2} \sum_{i,j} \sum_a \delta(A_{ij}^1, a) \log \nu_a^{\theta_i\theta_j}$ does not depend on the community structure. Similarly, by denoting $R_{ab}^{\theta_i\theta_j} = \log \frac{P_{ab}^{\theta_i\theta_j}}{Q_{ab}^{\theta_i\theta_j}}$, we find that

$$ 
\begin{aligned}
\log \mathbb{P}(A^t | A^{t-1}, \sigma, \theta) \\
&= \frac{1}{2} \sum_{i,j} \sum_{a,b} \delta(A_{ij}^{t-1}, a) \delta(A_{ij}^t, b) \left( \delta(\sigma_i, \sigma_j) R_{ab}^{\theta_i\theta_j} + \log Q_{ab}^{\theta_i\theta_j} \right) \\
&= \frac{1}{2} \sum_{i,j} \delta(\sigma_i, \sigma_j) \sum_{a,b} \delta(A_{ij}^{t-1}, a) \delta(A_{ij}^t, b) R_{ab}^{\theta_i\theta_j} + c_t(A),
\end{aligned}
$$

where $c_t(A) = \frac{1}{2} \sum_{i,j} \sum_{a,b} \delta(A_{ij}^{t-1}, a) \delta(A_{ij}^t, b) \log Q_{ab}^{\theta_i\theta_j}$ does not depend on the community structure. Simple

computations show that

$$ \sum_a \delta(A_{ij}^1, a) \rho_a^{\theta_i\theta_j} = A_{ij}^1 (\rho_1^{\theta_i\theta_j} - \rho_0^{\theta_i\theta_j}) + \rho_0^{\theta_i\theta_j} $$

and

$$ 
\begin{aligned}
\sum_{a,b} \delta(A_{ij}^{t-1}, a) \delta(A_{ij}^t, b) R_{ab}^{\theta_i\theta_j} \\
&= R_{00}^{\theta_i\theta_j} + A_{ij}^{t-1} (R_{10}^{\theta_i\theta_j} - R_{00}^{\theta_i\theta_j}) + A_{ij}^t (R_{01}^{\theta_i\theta_j} - R_{00}^{\theta_i\theta_j}) \\
&\quad + A_{ij}^{t-1} A_{ij}^t (R_{11}^{\theta_i\theta_j} - R_{01}^{\theta_i\theta_j} - R_{10}^{\theta_i\theta_j} + R_{00}^{\theta_i\theta_j}) \\
&= R_{00}^{\theta_i\theta_j} + A_{ij}^{t-1} \ell_{10}^{\theta_i\theta_j} + A_{ij}^t \ell_{01}^{\theta_i\theta_j} \\
&\quad + A_{ij}^{t-1} A_{ij}^t (\ell_{11}^{\theta_i\theta_j} - \ell_{01}^{\theta_i\theta_j} - \ell_{10}^{\theta_i\theta_j}).
\end{aligned}
$$

By collecting the above observations, we now find that $\log \mathbb{P}(A | \sigma, \theta)$ equals

$$ 
\begin{aligned}
\frac{1}{2} \sum_{\substack{i,j \\ \sigma_i = \sigma_j}} \left\{ A_{ij}^1 (\rho_1^{\theta_i\theta_j} - \rho_0^{\theta_i\theta_j}) + \rho_0^{\theta_i\theta_j} + (A_{ij}^1 - A_{ij}^T) \ell_{10}^{\theta_i\theta_j} \right\} \\
+ \frac{1}{2} \sum_{\substack{i,j \\ \sigma_i = \sigma_j}} \sum_{t=2}^{T} \left\{ (\ell_{01}^{\theta_i\theta_j} + \ell_{10}^{\theta_i\theta_j}) (A_{ij}^t - A_{ij}^{t-1} A_{ij}^t) \right. \\
\left. + \ell_{11}^{\theta_i\theta_j} A_{ij}^{t-1} A_{ij}^t - \log \frac{Q_{00}^{\theta_i\theta_j}}{P_{00}^{\theta_i\theta_j}} \right\} + c(A),
\end{aligned}
$$

where $c(A) = \sum_t c_t(A)$ does not depend on $\sigma$. Hence the claim follows. $\square$

### B. Sparse MLE and Modularity

Recall the structural assumptions (3)–(4) about the degree correction parameters. Because $P_{01}, Q_{01} = o(1)$, a first-order Taylor expansion yields

$$ \log \frac{1 - \theta_i \theta_j Q_{01}}{1 - \theta_i \theta_j P_{01}} = \theta_i \theta_j (P_{01} - Q_{01}) + o(P_{01} + Q_{01}), $$

as well as $\ell_{01}^{\theta_i\theta_j} \approx \log \frac{P_{01}}{Q_{01}}$, $\ell_{10}^{\theta_i\theta_j} \approx \log \frac{1-P_{11}}{1-Q_{11}}$ and $\ell_{11}^{\theta_i\theta_j} \approx \log \frac{P_{11}}{Q_{11}}$. Using these approximations in the MLE expression leads to the maximisation of

$$ \sum_{t=2}^{T} \sum_{i,j:\sigma_i=\sigma_j} (\tilde{a}_{ij}^t - \theta_i \theta_j (P_{01} - Q_{01})), \quad (19) $$

where $\tilde{a}_{ij}^t = \alpha(A_{ij}^t) + \beta(A_{ij}^t)$. Since $\mu$ and $\nu$ are stationary distributions,

$$ \mathbb{E} (A_{\text{new}}^t)_{ij} = \begin{cases} \theta_i \theta_j \mu_1 (1 - P_{11}) & \text{if } \sigma_i = \sigma_j, \\ \theta_i \theta_j \nu_1 (1 - Q_{11}) & \text{otherwise,} \end{cases} $$

$$ \mathbb{E} (A_{\text{pers}}^t)_{ij} = \begin{cases} \theta_i \theta_j \mu_1 P_{11} & \text{if } \sigma_i = \sigma_j, \\ \theta_i \theta_j \nu_1 Q_{11} & \text{otherwise.} \end{cases} $$
Therefore, using $W_{ij} = \sum_{t=2}^{T} \tilde{a}_{ij}$, we have

$$ \mathbb{E}W_{ij} = \begin{cases} (T-1)\theta_i\theta_j\mu_1(\alpha(1-P_{11}) + \beta P_{11}) & \text{if } \sigma_i = \sigma_j \\ (T-1)\theta_i\theta_j\nu_1(\alpha(1-Q_{11}) + \beta Q_{11}) & \text{otherwise.} \end{cases} $$

Since the community labelling is sampled uniformly at random, using the normalization for the $\theta_i$'s, we obtain $\bar{d}_i = (T-1)\theta_i N \frac{\mu_1(\alpha(1-P_{11})+\beta P_{11})+(K-1)\nu_1(\alpha(1-Q_{11})+\beta Q_{11})}{K}$,

together with $\bar{m} =$

$$ \frac{(T-1)N^2 \mu_1(\alpha(1-P_{11})+\beta P_{11})+(K-1)\nu_1(\alpha(1-Q_{11})+\beta Q_{11})}{2K} $$

we observe that

$$ \theta_i \theta_j (P_{01} - Q_{01}) = \gamma \frac{\bar{d}_i \bar{d}_j}{2\bar{m}}, $$

where $\gamma = (P_{01} - Q_{01})(T-1) \frac{\mu_1(\alpha(1-P_{11})+\beta P_{11})+(K-1)\nu_1(\alpha(1-Q_{11})+\beta Q_{11})}{K}$. We end the proof by employing (19). $\square$

### C. Modularity and Spectral Clustering

The regularized modularity of a partition $\sigma \in [K]^N$ of the graph A is defined as

$$ \mathcal{M}(A, \sigma, \gamma) = \sum_{i,j} \delta(\sigma_i, \sigma_j) \left( A_{ij} - \gamma \frac{d_i d_j}{2m} \right), $$

where $d = A1_n$ is a vector of degrees and $\gamma$ is a resolution parameter. The above equation can be rewritten as

$$ \mathcal{M}(A, \sigma, \gamma) = \operatorname{Tr} \tilde{Z}^T \left( A - \gamma \frac{dd^T}{2m} \right) \tilde{Z}, $$

where $\tilde{Z} \in \{0, 1\}^{N \times K}$ is the membership matrix associated with the vector $\sigma$, that is $\tilde{Z}_{ik} = 1$ for $k = \sigma_i$, and $\tilde{Z}_{ik} = 0$ otherwise. As maximising the modularity over $\sigma \in [K]^N$ is in general NP-complete [23], it is convenient to perform a continuous relaxation. Following [24], we transform the problem into

$$ \hat{X} = \arg \max_{\substack{X \in \mathbb{R}^{N \times K} \\ X^T DX = I_K}} \mathrm{Tr} X^T \left( A - \gamma \frac{dd^T}{2m} \right) X. \quad (20) $$

The predicted membership matrix $\hat{Z}$ is then recovered by performing an approximated solution to the following k-means problem (see [38])

$$ (\hat{Z}, \hat{Y}) = \arg \min_{Z \in \mathcal{Z}_{N,K}, Y \in \mathbb{R}^{K \times K}} \|ZY - \hat{X}\|_F. \quad (21) $$

The Lagrangian associated to the optimisation problem (20) is

$$ \mathrm{Tr} X^T \left( A - \gamma \frac{dd^T}{2m} \right) X - \mathrm{Tr} \left( \Lambda^T (X^T X - I_K) \right), $$

where $\Lambda \in \mathbb{R}^{K \times K}$ is a symmetric matrix of Lagrangian multipliers. Up to a change of basis, we can assume that $\Lambda$ is diagonal. The stationarity condition for the solution of (20) becomes

$$ (A - \gamma \frac{dd^T}{2m}) X = DX\Lambda \quad \text{and} \quad X^T DX = I_K, $$

which is a generalized eigenvalue problem: the columns of X are the generalized eigenvectors, and the diagonal elements of

$\Lambda$ are the eigenvalues. In particular, since the constant vector $1_n$ satisfies $(A - \gamma \frac{dd^T}{2m}) 1_n = (1 - \gamma) D1_n$, we conclude that the eigenvalues should be larger than $1 - \gamma$ for the partition to be meaningful.

Multiplying the first equation by $1_n^T$ leads to $(1 - \gamma) d^T X = d^T X \Lambda$, and therefore $d^T X = 0$ (using the previous remark on $\Lambda$). The system then simplifies in

$$ AX = DX\Lambda \quad \text{and} \quad X^T DX = I_K. $$

Defining a re-scaled vector $U = D^{-1/2}X$ shows that $U$ satisfies $D^{-1/2}AD^{-1/2}U = U\Lambda$ and $U^T U = I_K$. Thus, the columns of $U$ are eigenvectors of $D^{-1/2}AD^{-1/2}$ associated with the $K$ largest eigenvalues (or equivalently, the eigenvectors of $\mathcal{L} = I_N - D^{-1/2}AD^{-1/2}$ associated to the $K$ smallest eigenvalues).

## APPENDIX B

### PROOFS FOR SECTION III

**Proof of Proposition 3.1:** By letting $\Gamma \in \mathbb{R}^{K \times K}$ be the diagonal matrix whose diagonal entries $\Gamma_{kk}$ are the Lagrange multipliers associated with the constraints $(X^TDX)_{kk} = 1$, the Lagrangian of the relaxed optimisation problem (14) is

$$ \mathrm{Tr} (X^T W_{\tau} X + 2\lambda S^T X) + \mathrm{Tr} ((X^T DX - I_k) \Gamma). $$

This leads to the constrained linear system

$$ \begin{cases} W_{\tau}X + DX\Gamma = -\lambda S, \\ (X^TDX)_{kk} = 1, \end{cases} \quad (22) $$

whose unknowns are $\Gamma$ and $X$.

Firstly, we note that if $(\Gamma^1, X^1)$ and $(\Gamma^2, X^2)$ are solutions of the system (22), then (see Lemma B.1 for the computations)

$$ \begin{aligned} & C(X^1) - C(X^2) \\ &= \mathrm{Tr} \left[ (\Gamma^1 - \Gamma^2)(X^1 - X^2)^T D(X^1 - X^2) \right] \\ &= \sum_{k \in [K]} (\Gamma_{kk}^2 - \Gamma_{kk}^1) \|X_k^1 - X_k^2\|_D \end{aligned} $$

where $C(x) = \mathrm{Tr}(X^T W_{\tau} X + 2\lambda S^T X)$ is the cost function maximised in (14) and $\|x\|_D = x^T D x$. Hence, among the solution pairs $(\Gamma, X)$ of the system (22), the solution of the minimisation problem (14) is the vector $X$ associated with the diagonal matrix $\Gamma$ whose entries are the smallest.

Secondly, the eigenvalue decomposition of $-D^{-1/2}W_{\tau} D^{-1/2}$ reads as

$$ -D^{-1/2}W_{\tau}D^{-1/2} = Q\Delta Q^T, $$

where $\Delta = \mathrm{diag}(\delta_1, \dots, \delta_N)$ with $\delta_1 \le \dots \le \delta_N$ and $Q^T Q = I_N$. Therefore, after the change of variables $U = Q^T D^{1/2} X$ and $B = \lambda Q^T D^{-1/2} S$, the system (22) is transformed to

$$ \begin{cases} \Delta U = U\Gamma + B, \\ (U^T U)_{kk} = 1 \quad \forall k \in [K]. \end{cases} $$

Thus, the solution $\hat{X} = (\hat{X}_{.1}, \dots, \hat{X}_{.K})$ of the optimisation problem (14) satisfies

$$ (-W_{\tau} - \gamma_k^* D) \hat{X}_{.k} = \lambda S_{.k}, $$
where for all $k \in [K]$, $\gamma_k^*$ is the smallest solution of the explicit secular equation [39]

$$ \sum_{i=1}^{N} \left( \frac{B_{ik}}{\delta_i - \gamma} \right)^2 - 1 = 0. \quad \square $$

Lemma B.1: If $(\Gamma^1, X^1)$ and $(\Gamma^2, X^2)$ are solutions of the system (22), then

$$ C(X^1) - C(X^2) = \sum_{k \in [K]} (\Gamma_{kk}^2 - \Gamma_{kk}^1) \|X_{\cdot k}^1 - X_{\cdot k}^2\|_D, $$

where $C(X) = \text{Tr}(X^T W_{\tau} X + 2\lambda S^T X)$ and $\|x\|_D = x^T D x$.

Proof: Because $(\Gamma^1, X^1)$ and $(\Gamma^2, X^2)$ are solutions of (22), it holds that

$$ \begin{aligned} -W_{\tau}X^1 - DX^1\Gamma^1 &= \lambda S, \\ -W_{\tau}X^2 - DX^2\Gamma^2 &= \lambda S. \end{aligned} \quad (23) $$

Using the definition of $C(X)$ and the fact that $(\Gamma^1, X^1)$ is solution of (22), we notice that

$$ C(X^1) = \text{Tr} \left( -(X^1)^T DX^1 \Gamma^1 + \lambda (X^1)^T S \right), $$

and similarly for $C(X_2)$. Moreover, multiplying the first equation of (23) (resp., the second equation) by $(X^2)^T$ (resp., by $(X^1)^T$) gives

$$ \begin{aligned} \lambda(X^2)^T S &= -(X^2)^T W_{\tau} X^1 - (X^2)^T DX^1 \Gamma^1, \\ \lambda(X^1)^T S &= -(X^1)^T W_{\tau} X^2 - (X^1)^T DX^2 \Gamma^2. \end{aligned} $$

Because $(X^1)^T W_{\tau} X^2 = (X^2)^T W_{\tau} X^1$, we have

$$ \begin{aligned} \text{Tr} (\lambda S^T (X^1 - X^2)) &= \text{Tr}(-(X^1)^T DX^2 \Gamma^2 + (X^2)^T DX^1 \Gamma^1). \end{aligned} $$

Therefore,

$$ C(X^1) - C(X^2) = \text{Tr} \left( -(X^1)^T DX^1 \Gamma^1 + (X_2)^T DX^1 \Gamma^1 \right) \\ - \text{Tr} \left( -(X^2)^T DX^2 \Gamma^2 + (X^1)^T DX^2 \Gamma^2 \right). \quad (24) $$

Moreover, because $\Gamma^1$ and $D$ are symmetric, we have

$$ \begin{aligned} \text{Tr} \left( -(X^1)^T DX^1 \Gamma^1 + (X_2)^T DX^1 \Gamma^1 \right) &= \sum_{k \in [K]} \Gamma_{kk}^1 \left( \left( -(X^1)^T DX^1 \right)_{kk} + \left( (X_2)^T DX^1 \right)_{kk} \right) \\ &= \sum_{k \in [K]} \Gamma_{kk}^1 \left( -(X_{\cdot k}^1)^T DX_{\cdot k}^1 + (X_{\cdot k}^2)^T DX_{\cdot k}^1 \right). \end{aligned} $$

The constraints $X_{\cdot k}^T DX_{\cdot k} = 1$ verified by both $X_1$ and $X_2$ leads to

$$ \begin{aligned} & -(X_{\cdot k}^1)^T DX_{\cdot k}^1 + (X_{\cdot k}^2)^T DX_{\cdot k}^1 \\ &= - \frac{(X_{\cdot k}^1)^T DX_{\cdot k}^1 + (X_{\cdot k}^2)^T DX_{\cdot k}^2}{2} + (X_{\cdot k}^2)^T DX_{\cdot k}^1 \\ &= - \|X_{\cdot k}^1 - X_{\cdot k}^2\|_D, \end{aligned} $$

where the last equality uses $(X_{\cdot k}^2)^T DX_{\cdot k}^1 = (X_{\cdot k}^1)^T DX_{\cdot k}^2$. This ensures that

$$ \begin{aligned} \text{Tr} \left( -(X^1)^T DX^1 \Gamma^1 + (X_2)^T DX^1 \Gamma^1 \right) &= - \sum_{k \in [K]} \Gamma_{kk}^1 \|X_{\cdot k}^1 - X_{\cdot k}^2\|_D. \end{aligned} \quad (25) $$

Similarly, we can establish that

$$ \begin{aligned} \text{Tr} \left( -(X^2)^T DX^2 \Gamma^2 + (X^1)^T DX^2 \Gamma^2 \right) &= - \sum_{k \in [K]} \Gamma_{kk}^2 \|X_{\cdot k}^1 - X_{\cdot k}^2\|_D. \end{aligned} \quad (26) $$

We finish the proof by combining (25) and (26) with (24). □

REFERENCES

[1] S. Fortunato and M. E. Newman, “20 years of network community detection,” *Nature Phys.*, vol. 18, no. 8, pp. 848–850, 2022.

[2] P. Holme and J. Saramäki, “Temporal networks,” *Phys. Rep.*, vol. 519, no. 3, pp. 97–125, 2012.

[3] S. Paul and Y. Chen, “Consistent community detection in multi-relational data through restricted multi-layer stochastic blockmodel,” *Electron. J. Statist.*, vol. 10, no. 2, pp. 3807–3870, 2016.

[4] K. Avrachenkov, M. Dreveton, and L. Leskelä, “Community recovery in non-binary and temporal stochastic block models,” 2022, *arXiv:2008.04790v5*.

[5] S. Bhattacharyya and S. Chatterjee, “General community detection with optimal recovery conditions for multi-relational sparse networks with dependent layers,” 2020, *arXiv:2004.03480*.

[6] J. Lei and K. Z. Lin, “Bias-adjusted spectral clustering in multi-layer stochastic block models,” *J. Amer. Stat. Assoc.*, pp. 1–13, 2022.

[7] A. Ghasemian, P. Zhang, A. Clauset, C. Moore, and L. Peel, “Detectability thresholds and optimal algorithms for community structure in dynamic networks,” *Phys. Rev. X*, vol. 6, no. 3, 2016, Art. no. 031005.

[8] C. Matias and V. Miele, “Statistical clustering of temporal networks through a dynamic stochastic block model,” *J. Roy. Stat. Soc.: Ser. B.* (*Statistical Methodol.*), vol. 79, no. 4, pp. 1119–1141, 2017.

[9] Y. Chi, X. Song, D. Zhou, K. Hino, and B. L. Tseng, “Evolutionary spectral clustering by incorporating temporal smoothness,” in *Proc. 13th ACM SIGKDD Int. Conf. Knowl. Discov. Data Mining*, 2007, pp. 153–162.

[10] F. Liu, D. Choi, L. Xie, and K. Roeder, “Global spectral clustering in dynamic networks,” *Proc. Nat. Acad. Sci.*, vol. 115, no. 5, pp. 927–932, 2018.

[11] L. Dall'Amico, R. Couillet, and N. Tremblay, “Community detection in sparse time-evolving graphs with a dynamical bethe-hessian,” in *Proc. Adv. Neural Inf. Process. Syst.*, 2020, pp. 7486–7497.

[12] 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.

[13] A. R. Pamfil, S. D. Howison, R. Lambiotte, and M. A. Porter, “Relating modularity maximization and stochastic block models in multilayer networks,” *SIAM J. Math. Data Sci.*, vol. 1, no. 4, pp. 667–698, 2019.

[14] G. Rossetti and R. Cazabet, “Community discovery in dynamic networks: A survey,” *ACM Comput. Surv.*, vol. 51, no. 2, pp. 1–37, Feb. 2018.

[15] K. Avrachenkov and M. Dreveton, *Statistical Analysis of Networks*. Boston, MA, USA: Now, 2022.

[16] P. Barucca, F. Lillo, P. Mazzarisi, and D. Tantari, “Disentangling group and link persistence in dynamic stochastic block models,” *J. Stat. Mechanics: Theory Exp.*, vol. 2018, 2018, Art. no. 123407.

[17] P. Holland, K. B. Laskey, and S. Leinhardt, “Stochastic blockmodels: First steps,” *Social Netw.*, vol. 5, pp. 109–137, 1983.

[18] B. Karrer and M. E. J. Newman, “Stochastic blockmodels and community structure in networks,” *Phys. Rev. E*, vol. 83, Jan. 2011, Art. no. 016107.

[19] K. Avrachenkov, M. Dreveton, and L. Leskelä, “Recovering communities in temporal networks using persistent edges,” in *Proc. Comput. Data Social Netw.*, *CSoNet*, 2021, pp. 243–254.

[20] M. E. Newman and M. Girvan, “Finding and evaluating community structure in networks,” *Phys. Rev. E*, vol. 69, no. 2, 2004, Art. no. 026113.

[21] J. Reichardt and S. Bornholdt, “Statistical mechanics of community detection,” *Phys. Rev. E*, vol. 74, no. 1, 2006, Art. no. 016110.
[22] M. E. Newman, "Equivalence between modularity optimization and maximum likelihood methods for community detection," *Phys. Rev. E*, vol. 94, no. 5, 2016, Art. no. 052315.

[23] U. Brandes et al., "On finding graph clusterings with maximum modularity," in *Proc. Int. Workshop Graph-Theoretic Concepts Comput. Sci.*, 2007, pp. 121-132.

[24] M. E. J. Newman, "Spectral methods for community detection and graph partitioning," *Phys. Rev. E*, vol. 88, Oct. 2013, Art. no. 042822.

[25] M. Rosvall, A. V. Esquivel, A. Lancichinetti, J. D. West, and R. Lambiotte, "Memory in network flows and its effects on spreading dynamics and community detection," *Nature Commun.*, vol. 5, no. 1, 2014, Art. no. 4630.

[26] J. Xu, T. L. Wickramarathne, and N. V. Chawla, "Representing higher-order dependencies in networks," *Sci. Adv.*, vol. 2, no. 5, 2016, Art. no. e1600028.

[27] L. Eldén and H. Park, "A procrustes problem on the stiefel manifold," *Numerische Mathematik*, vol. 82, no. 4, pp. 599-619, 1999.

[28] L.-H. Zhang, W. H. Yang, C. Shen, and J. Ying, "An eigenvalue-based method for the unbalanced procrustes problem," *SIAM J. Matrix Anal. Appl.*, vol. 41, no. 3, pp. 957-983, 2020.

[29] T. Yang, Y. Chi, S. Zhu, Y. Gong, and R. Jin, "Detecting communities and their evolutions in dynamic social networks-a Bayesian approach," *Mach. Learn.*, vol. 82, no. 2, pp. 157-189, 2011, doi: 10.1007/s10994-010-5214-7.

[30] K. S. Xu and A. O. Hero, "Dynamic stochastic blockmodels for time-evolving social networks," *IEEE J. Sel. Topics Signal Process.*, vol. 8, no. 4, pp. 552-562, Aug. 2014.

[31] M. Pensky and T. Zhang, "Spectral clustering in the dynamic stochastic block model," *Electron. J. Statist.*, vol. 13, no. 1, pp. 678-709, 2019, doi: 10.1214/19-EJS1533.

[32] K. Xu, "Stochastic block transition models for dynamic networks," in *Proc. Int. Conf. Artif. Intell. Statist.*, 2015, pp. 1079-1087. [Online]. Available: https://proceedings.mlr.press/v38/xu15.html

[33] J. Fournet and A. Barrat, "Contact patterns among high school students," *PLoS One*, vol. 9, no. 9, pp. 1-17, 2014.

[34] R. Mastrandrea, J. Fournet, and A. Barrat, "Contact patterns in a high school: A comparison between data collected using wearable sensors, contact diaries and friendship surveys," *PLoS One*, vol. 10, no. 9, pp. 1-26, 2015.

[35] J. Stehlé et al., "High-resolution measurements of face-to-face contact patterns in a primary school," *PLoS One*, vol. 6, no. 8, 2011, Art. no. e23176.

[36] M. Génois, C. L. Vestergaard, J. Fournet, A. Panisson, I. Bonmarin, and A. Barrat, "Data on face-to-face contacts in an office building suggest a low-cost vaccination strategy based on community linkers," *Netw. Sci.*, vol. 3, no. 3, pp. 326-347, 2015.

[37] M. Génois and A. Barrat, "Can co-location be used as a proxy for face-to-face contacts?," *EPJ Data Sci.*, vol. 7, no. 1, pp. 1-18, 2018.

[38] A. Kumar and R. Kannan, "Clustering with spectral norm and the k-means algorithm," in *Proc. IEEE 51st Annu. Symp. Foundations Comput. Sci.*, 2010, pp. 299-308.

[39] W. Gander, G. H. Golub, and U. Von Matt, "A constrained eigenvalue problem," *Linear Algebra Appl.*, vol. 114, pp. 815-839, 1989.

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

**Konstantin Avrachenkov** received the master's degree in control theory from St. Petersburg State Polytechnic University, Saint Petersburg, Russia, in 1996, the Ph.D. degree in mathematics from the University of South Australia, Adelaide, Australia, in 2000, and the Habilitation (Doctor of Science) degree from the University of Nice Sophia Antipolis, Nice, France, in 2010. He is currently the Director of Research with INRIA, Sophia Antipolis, France. His research interests include Markov chains, Markov decision processes, stochastic games and singular perturbations. He applies these methodological tools to the modelling and control of networks and to design data mining and machine learning algorithms. He was the recipient of the five best paper awards. He is an Associate Editor for *Probability in the Engineering and Informational Sciences*, ACM TOMPECS, Stochastic Models, *IEEE Network Magazine* and on the Advisory Board of the *International Journal of Performance Evaluation*.

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

**Maximilien Dreveton** received the bachelor's and master's degrees in physics from the Ecole Normale Supérieure de Lyon, Lyon, France in 2013 and 2015, respectively, and the Ph.D. degree in computer science from Inria Sophia Antipolis, Biot, France, in 2022. He is currently a Postdoctoral Fellow with École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland. His research interests include statistical analysis of random graphs, and more particularly community detection in networks.

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

**Lasse Leskelä** received the M.Sc. and D.Sc. degrees in mathematics from the Helsinki University of Technology, Espoo, Finland, in 1999 and 2005, respectively. He also holds a title of Docent with the University of Jyväskylä, Jyväskylä, Finland. He is currently an Associate Professor of applied mathematics with Aalto University, Espoo. His research interests include network statistics, random graphs, probability theory, and stochastic processes.