SAN 阅读笔记
目录

第 06 章精校翻译:时序网络社区检测

第 6 章 时序网络中的社区检测(Community Detection in Temporal Networks)

前面各章集中于静态相互作用的研究,相互作用由一个二值数或一个正权重表示。然而,在许多应用领域中,相互作用随时间变化。这类数据集的纵向(longitudinal)本质要求我们用时序网络模型——由张量表示的时序网络(temporal network)——取代经典的基于图的模型(Holme and Saramäki, 2012; Kivelä et al., 2014)。我们注意到,不谨慎地处理时间方面——例如沿时间轴对时序数据做聚合或平滑——可能导致宝贵信息的丢失。

时序网络中的社区检测(community detection)问题近年来吸引了科学界相当多的关注。它在产生许多有趣结果的同时,也导致了彼此歧异的术语与算法的爆炸式增长。

在本章中,我们首先把现有的带社区时序网络模型统一到一个框架中。然后,我们将说明现有工作可以分为两大类:社区成员固定的模型与社区成员随时间变化的模型。随后我们分别研究这两种情形。

6.1 带社区的时序网络的一般模型(A General Model of Temporal Networks with Communities)

6.1.1 成员结构与相互作用结构(Membership and Interaction Structures)

我们考虑一个具有 $n$ 个节点、$K$ 个块和 $T$ 个时间快照(snapshot)的时序网络分块模型。观测数据由 $T$ 个邻接矩阵的列表 $(A^1, \ldots, A^T)$ 组成,其中每个矩阵 $A^t \in \{0,1\}^{n\times n}$ 描述网络在某一特定时刻的快照。此外,在时刻 $t$,节点集被划分为 $K$ 个潜在社区,我们用 $Z_{it}$ 表示节点 $i$ 在时刻 $t$ 的标签。

矩阵 $Z \in [K]^{n\times T}$ 表示成员结构(membership structure)。每一列 $Z_{\cdot t} \in [K]^n$ 由给定时刻 $t$ 各节点的社区标签组成,而行 $Z_{i\cdot} \in [K]^T$ 表示节点 $i$ 的成员模式(membership pattern,即节点 $i$ 的社区标签的演化)。

我们假设各节点的成员模式相互独立,且服从 $[K]^T$ 上的一个概率分布 $\mathscr{P}$。因此,

$$ \mathbb{P}(Z) = \prod_{i=1}^{n} \mathscr{P}\left(Z_{i\cdot}\right).\tag{6.1} $$

在块成员结构 $Z$ 的条件下,我们要生成一个随机张量 $A = \left(A_{ij}^{t}\right) \in \{0,1\}^{n\times n\times T}$,它以节点对 $\{i,j\}$ 与时刻 $t$ 为指标,满足 $A_{ij}^{t} = A_{ji}^{t}$,并使得节点对之间的模式相互作用相互独立。我们用 $B_{k^{1:T},\ell^{1:T}}\left(x^{1:T}\right)$ 表示在一对分别具有块模式 $k^{1:T} = (k_1, \ldots, k_T) \in [K]^T$ 与 $\ell^{1:T} = (\ell_1, \ldots, \ell_T) \in [K]^T$ 的节点之间,观测到相互作用模式 $x^{1:T} \in \{0,1\}^T$ 的概率。相互作用结构(interaction structure) $B$ 因而是概率测度的一个族 $B = \left(B_{k^{1:T},\ell^{1:T}}\right)$。在 $B$ 与 $Z$ 的条件下,随机张量 $A$ 的分布定义为

$$ \mathbb{P}(A \mid Z, B) = \prod_{1 \leq i < j \leq n} B_{Z_{i\cdot}, Z_{j\cdot}}\left(A_{ij}^{1:T}\right).\tag{6.2} $$

模型 (6.1)–(6.2) 是具有 $n$ 个节点、$K$ 个簇、$T$ 个快照的时序网络分块模型最一般的表达。由于块成员结构的大小为 $K \times T$,概率测度 $B_{k^{1:T},\ell^{1:T}}$ 共有 $\frac{(KT)^2}{2}$ 种选择。保持这种完全的普遍性会导致一个过于复杂的模型。下一节详述若干有趣的特例。

6.1.2 时序网络模型的例子(Examples of Temporal Network Models)

静态成员、动态相互作用(Static memberships, dynamic interactions) 若矩阵 $Z$ 的各列相等,则称块成员结构 $Z$ 是静态的。等价地,每个节点的社区标记不随时间变化。此时我们可以简单地用一个向量 $z \in [K]^n$ 表示静态社区标记。进一步,块相互作用结构 $B$ 退化为一个相互作用核(interaction kernel)$f = (f_{k\ell})_{k,\ell\in[K]}$,它是 $\mathcal{S} = \{0,1\}^T$ 上的一族概率分布,满足 $f_{k\ell} = f_{\ell k}$。这定义了对称相互作用张量 $A \in \mathcal{S}^{n\times n}$ 的概率分布

$$ \mathbb{P}(A \mid z) = \prod_{1 \leq i < j \leq n} f_{z_i z_j}\left(A_{ij}\right)\tag{6.3} $$

它表示一个具有静态社区结构的时序分块模型。若相互作用核形如

$$ f_{k\ell} = \left\{ \begin{array}{ll} f_{\mathrm{in}}, & \quad \text{若 } k = \ell, \\ f_{\mathrm{out}}, & \quad \text{其他}, \end{array} \right. $$

则称模型是同质(homogeneous)的。这里 $f_{\mathrm{in}}$ 表示块内相互作用的分布,而跨块的相互作用服从 $f_{\mathrm{out}}$。

例 6.1 时间独立的相互作用

设 $x = (x_1, \ldots, x_T) \in \{0,1\}^T$。若对所有 $k, \ell \in [K]$,有 $f_{k\ell}(x) = \prod_{t=1}^{T} \mu_{k\ell}(x_t)$,其中 $\mu_{k\ell}$ 是 $\{0,1\}$ 上的概率分布,则称具有静态成员结构的时序 SBM 具有时间独立的相互作用(temporally independent interactions)。具有静态成员结构与时间独立相互作用的时序分块模型,对应于一个二值 SBM 的 $T$ 次独立观测。

例 6.2 马尔可夫相互作用

设 $x = (x_1, \ldots, x_T) \in \{0,1\}^T$。若对所有 $k, \ell \in [K]$,有 $f_{k\ell} = \mu_{k\ell}(x_1) \prod_{t=2}^{T} P_{k\ell}(x_{t-1}, x_t)$,其中 $\mu_{k\ell}$ 是 $\{0,1\}$ 上的概率分布,$P_{k\ell}$ 是表示相邻两个快照之间转移概率的 $2 \times 2$ 随机矩阵,则称具有静态成员的时序 SBM 具有马尔可夫相互作用(Markov interactions)。若对所有 $k, \ell \in [K]$ 与 $a, b \in \{0,1\}$ 有 $P_{k\ell}(a, b) = \mu_{k\ell}(b)$,则回到例 6.1 所描述的模型。

时间独立的相互作用(Temporally independent interactions) 若在任意时刻 $t$,节点 $i$ 与 $j$ 之间的二值相互作用都按 $i$ 与 $j$ 在时刻 $t$ 的社区标记重新抽样,则称块相互作用结构是时间独立的。此时随机张量 $A$ 的分布由下式给出:

$$ \mathbb{P}(A \mid Z, B) = \prod_{1 \leq i < j \leq n} Q_{Z_{it}, Z_{jt}}\left(A_{ij}^{t}\right)\tag{6.4} $$

其中 $Q = (Q_{k\ell})_{k,\ell\in[K]}$ 是 $\{0,1\}$ 上的一族分布。

马尔可夫成员结构(Markov membership structure) 设 $z^{1:T} = (z_1, \dots, z_T) \in [K]^T$ 表示一个成员模式,并回忆 $\mathscr{P}$ 是节点社区指派的分布(见式 (6.1))。若

$$ \mathscr{P}(z) = \alpha_{z_1} \prod_{t=2}^{T} \pi_{z_{t-1}, z_t}, $$

其中 $\alpha$ 是 $[K]$ 上的初始概率分布、$\pi$ 是 $K \times K$ 转移概率矩阵,则称模型具有马尔可夫成员结构。

例 6.3 马尔可夫成员结构的一个参数化

设 $\mathbb{1}_K = (1, \ldots, 1)^T$ 为全 $1$ 的 $K \times 1$ 向量。由 $\alpha = \frac{1}{K}\mathbb{1}_K$ 与 $\pi = rI_K + \frac{1-r}{K}\mathbb{1}_K \mathbb{1}_K^T$ 定义的具有马尔可夫成员结构的模型,对应于如下模型:

  • 在初始时刻 $t = 1$,所有节点的社区标记独立地均匀随机选取;
  • 在时刻 $t \geq 2$,给定节点 $i$ 以概率 $r$ 留在它在时刻 $t-1$ 所在的社区,而以概率 $1-r$ 被指派到一个均匀随机选取的新社区。

6.2 静态社区成员的网络(Networks with Static Community Memberships)

6.2.1 马尔可夫相互作用 SBM 的恢复阈值(Recovery Thresholds in SBM with Markov Interaction)

若社区成员是静态的且时间相互作用是马尔可夫的(另见例 6.2),该模型称为马尔可夫随机分块模型(Markov Stochastic Block Model)。在同质马尔可夫 SBM 中,相互作用核形如

$$ \begin{array}{rl} & f_{\mathrm{in}} = \mu_{x_1} P_{x_1, x_2} \cdots P_{x_{T-1}, x_T}, \\ & f_{\mathrm{out}} = \nu_{x_1} Q_{x_1, x_2} \cdots Q_{x_{T-1}, x_T}, \end{array}\tag{6.5} $$

其中 $\mu, \nu$ 是 $\{0,1\}$ 上的初始概率分布,$P, Q$ 是 $\{0,1\}$ 上的转移概率矩阵。

在稀疏情形(sparse regime)下,任意特定节点对之间观测到非零相互作用的概率很小,即

$$ \max\left\{ \mu_1, \nu_1, P_{01}, Q_{01} \right\} \leq \rho.\tag{6.6} $$

一个特例是假设对某些常数 $u, v, p_{01}, q_{01} \in (0, \infty)$,有

$$ \mu_1 = u\rho, \quad \nu_1 = v\rho, \quad P_{01} = p_{01}\rho, \quad Q_{01} = q_{01}\rho.\tag{6.7} $$

在这一假设下,一个 $f$ 分布信号中 $1$ 的期望个数为 $\mathbb{E}\sum_{t=1}^{T} X_t \leq \mu_1 + (T-1)P_{01} = O(\rho T)$。因此,当 $\rho T = o(1)$ 时,在任意特定节点对中观测到相互作用的概率很小。

下述命题给出了当 $n \gg 1$ 且 $T \gg 1$ 时稀疏马尔可夫 SBM 的恢复条件,其证明可见 (Avrachenkov et al., 2022)。一致估计量与强一致估计量的概念已在 4.4.3 节定义。

命题 6.1 稀疏马尔可夫 SBM 的恢复条件

考虑一个由 $n \gg 1$ 个节点、$K \asymp 1$ 个块、$T \gg 1$ 个快照组成的同质马尔可夫 SBM,其中 $f_{\mathrm{in}}$ 与 $f_{\mathrm{out}}$ 是由 (6.5) 定义并满足 (6.7) 的马尔可夫链分布,稀疏参数 $\rho$ 满足 $\rho T \ll 1$。假设 $P_{11}$ 与 $Q_{11}$ 为常数,满足 $(P_{11}, Q_{11}) \neq (1, 1)$,且 $(p_{01}, P_{11}) \neq (q_{01}, Q_{11})$。令

$$ \tilde{I} = \left( \sqrt{p_{01}} - \sqrt{q_{01}} \right)^{2} + 2\sqrt{p_{01} q_{01}}\, H_{11}^{2}, $$

其中 $H_{11}^{2} = 1 - \dfrac{\sqrt{(1 - P_{11})(1 - Q_{11})}}{1 - \sqrt{P_{11} Q_{11}}}$ 是参数分别为 $P_{11}$ 与 $Q_{11}$ 的两个几何分布之间的平方 Hellinger 散度。则:

(i) 当 $\rho T \lesssim \frac{1}{n}$ 时一致估计量不存在,当 $\rho T \gg \frac{1}{n}$ 时一致估计量存在;

(ii) 当 $\rho T \ll \frac{\log n}{n}$ 时强一致估计量不存在,当 $\rho T \gg \frac{\log n}{n}$ 时强一致估计量存在;

(iii) 在临界情形(critical regime)$\rho T = (1 + o(1))\tau\frac{\log n}{n}$ 中($\tau$ 为某常数),当 $\tau\tilde{I} < K$ 时强一致估计量不存在,当 $\tau\tilde{I} > K$ 时强一致估计量存在。

Tips:这是本章与第 4 章的承接点:静态成员 + 马尔可夫交互把第 4 章 SBM 的恢复阈值中的信号强度从 $\rho$ 提升为 $\rho T$——快照数 $T$ 相对静态 SBM 带来信息增益(另见注 6.1、6.2)。原书未给证明,证明外包至 (Avrachenkov et al., 2022)。
查看学习笔记对该命题(书内无证明,证明见 Avrachenkov et al., 2022)的说明与条件梳理

量 $\rho T \tilde{I}$ 对应于两个马尔可夫链分布 $f_{\mathrm{in}}$ 与 $f_{\mathrm{out}}$ 之间 Rényi 散度的 Taylor 展开的主项。更多细节与证明见 (Avrachenkov et al., 2022)。

注 6.1 极稀疏情形下的恢复

回忆例 4.1:二值 SBM 中的一致恢复要求 $\rho \gg n^{-1}$。特别地,命题 6.1 表明,只要快照数足够大,即使在非常稀疏的情形下一致恢复也是可能的。例如,若 $\rho = \frac{1}{n}$,则 $T$ 至少需为 $\omega(1)$ 量级,一致恢复才有可能。

注 6.2 与静态 SBM 临界情形的对照

信号强度为 $\rho T = (1 + o(1))\tau\frac{\log n}{n}$ 的参数情形是一个有趣的临界情形。事实上,在该情形下,强一致的相变发生在 $\tau\left(\sqrt{p_{01}} - \sqrt{q_{01}}\right)^{2} + 2\tau\sqrt{p_{01}q_{01}}\,H_{11}^{2} > K$。作为对比,静态 SBM 中强一致的有趣参数情形是 $\rho = (1 + o(1))\frac{\log n}{n}$(见例 4.2)。

6.2.2 马尔可夫动力学的在线似然算法(Online Likelihood-based Algorithms for Markov Dynamics)

在本节中,我们为具有静态社区成员的时序网络聚类推导一个算法。我们考虑这样的情形:数据逐快照到达,并在每个时刻更新社区成员的在线估计。

模型参数已知(Model parameters are known) 给定 $A^{1:t} = \left(A^1, \cdots, A^t\right)$,定义对数似然比矩阵

$$ M_{ij}^{t} = \log \frac{f_{\mathrm{in}}\left(A_{ij}^{1:t}\right)}{f_{\mathrm{out}}\left(A_{ij}^{1:t}\right)},\tag{6.8} $$

其中 $f_{\mathrm{in}}$ 与 $f_{\mathrm{out}}$ 是块内与块间相互作用概率。特别地,给定节点标记 $z$ 时观测到图序列 $A^{1:t}$ 的概率的对数等于

$$ \log \mathbb{P}(A \mid z) = \frac{1}{2}\sum_{i}\sum_{j \neq i} M_{ij}^{t}\,\mathbb{1}(z_j = z_i) + \frac{1}{2}\sum_{i}\sum_{j \neq i} f_{\mathrm{out}}\left(A_{ij}^{1:t}\right). $$

因此,给定由前 $t-1$ 个快照的观测计算出的指派 $\hat{z}^{t-1}$,可以计算一个新的指派 $\hat{z}^{t}$,使得节点 $i$ 被指派到使下式最大的任意块 $k$:

$$ \mathcal{L}_{i,k}^{t} = \sum_{j \neq i} M_{ij}^{t}\,\delta_{\hat{z}_{j}^{t-1} k}.\tag{6.9} $$

只有当 $M^t$ 能容易地由 $M^{t-1}$ 算出时,这个公式才有意义。马尔可夫演化正是这种情形。事实上,若 $f_{\mathrm{in}}$ 与 $f_{\mathrm{out}}$ 由 (6.5) 给出,则式 (6.8) 定义的累积对数似然矩阵可以按 $M^t = M^{t-1} + \Delta^t$ 递归计算,其中

$$ M_{ij}^{1} = \log \frac{\mu}{\nu}\left(A_{ij}^{1}\right) \quad \text{且} \quad \Delta_{ij}^{t} = \log \frac{P}{Q}\left(A_{ij}^{t-1}, A_{ij}^{t}\right). $$

我们将此总结为算法 15。需要强调的是,该算法以在线自适应的方式工作。

算法 15 的时间复杂度(最坏情形复杂度)为 $O(Kn^2T)$,再加上初始聚类的时间复杂度。空间复杂度为 $O(n^2)$。

算法 15 块相互作用参数已知时同质马尔可夫动力学的在线聚类

输入:相互作用张量 $(A_{ij}^{t})$;块相互作用参数 $\mu, \nu, P, Q$;社区个数 $K$;静态图聚类算法(记为 $\mathsf{algo}$)。

输出:节点标记 $\hat{z} = \left(\hat{z}_1, \dots, \hat{z}_n\right) \in [K]^n$。

  1. 初始化:计算 $\hat{z} \gets \mathsf{algo}(A^1)$,并对 $i, j = 1, \dots, n$ 令 $M_{ij} \gets \log \frac{\mu}{\nu}\left(A_{ij}^{1}\right)$。
  2. for $t = 2, \ldots, T$ do
  3. 对 $i, j = 1, \dots, n$ 计算 $\Delta_{ij} \gets \log \frac{P}{Q}\left(A_{ij}^{t-1}, A_{ij}^{t}\right)$;
  4. 更新 $M \gets M + \Delta$;
  5. for $i = 1, \ldots, n$ do
  6. 对 $k = 1, \ldots, K$ 令 $L_{ik} \gets \sum_{j \neq i} M_{ij}\,\delta_{\hat{z}_j k}$;
  7. 令 $\hat{z}_i \gets \arg\max_{1 \le k \le K} L_{ik}$。

返回:$\hat{z}$。

此外,我们注意到:

  • 由于在每个时刻 $\Delta_{ij}$ 只能取四个值之一,这四个不同的 $\Delta_{ij}$ 值可以预先计算并存储,以避免计算 $n^2 T$ 次对数;
  • $n \times K$ 矩阵 $(L_{ik})$ 可以作为矩阵乘积 $L = M^0 Z$ 计算,其中 $M^0$ 是把 $M$ 的对角线置零所得的矩阵,$Z$ 是 $\hat{z}$ 的独热表示(one-hot representation),即当 $\hat{z}_i = k$ 时 $Z_{ik} = 1$,否则为 $0$;
  • 对于稀疏网络,通过忽略 $0 \to 0$ 转移并只存储非零元,时间与空间复杂度(平均复杂度)可以降低一个 $d/n$ 因子,其中 $d$ 是单个快照中的平均节点度。

参数未知时的推广(Extension when the parameters are unknown) 算法 15 要求先验地知道相互作用参数。实践中往往并非如此,必须在恢复社区的过程中学习参数。在本小节中,我们通过在运行中估计参数来改造算法 15。

设 $n_{ab}(i,j)$ 是节点 $i$ 与 $j$ 之间相互作用模式中观测到的 $a \to b$ 转移次数,并令 $n_a(i,j) = \sum_b n_{ab}(i,j)$。设 $P(i,j)$ 是节点对 $(i,j)$ 的模式相互作用演化的 $2 \times 2$ 转移概率矩阵。由大数定律(针对平稳且遍历的随机过程),经验转移概率

$$ \widehat{P}_{ab}(i,j) = \frac{n_{ab}(i,j)}{n_a(i,j)}\tag{6.10} $$

在 $T \gg 1$ 时以高概率接近 $P(i,j)$。

$P$ 的一个估计量可以通过在被预测为属于同一社区的节点对上对这些概率取平均而得到。更准确地说,在观测了 $t$ 个快照($t \geq 2$)之后,给定预测的社区指派 $\hat{z}^t$,对 $a, b \in \{0,1\}$ 定义

$$ \widehat{P}_{ab}^{t} = \frac{1}{\left| \left\{ (i,j) : \hat{z}_i^{t} = \hat{z}_j^{t} \right\} \right|} \sum_{(i,j) : \hat{z}_i^{t} = \hat{z}_j^{t}} \frac{n_{ab}^{t}(i,j)}{n_a^{t}(i,j)},\tag{6.11} $$

其中

$$ n_{ab}^{t}(i,j) = \sum_{t' = 1}^{t-1} \mathbb{1}\left(A_{ij}^{t'} = a\right)\mathbb{1}\left(A_{ij}^{t'+1} = b\right) $$

是在前 $t$ 个快照中见到的、节点 $i$ 与 $j$ 之间相互作用模式里 $a \to b$ 转移的次数($a, b \in \{0,1\}$),而

$$ n_a^{t}(i,j) = \sum_{b=0}^{1} n_{ab}^{t}(i,j). $$

类似地,

$$ \widehat{Q}_{ab}^{(t)} = \frac{1}{\left| \left\{ (i,j) : \hat{z}_i^{t} \neq \hat{z}_j^{t} \right\} \right|} \sum_{(i,j) : \hat{z}_i^{t} \neq \hat{z}_j^{t}} \frac{n_{ab}^{(t)}(i,j)}{n_a^{(t)}(i,j)},\tag{6.12} $$

是 $Q_{ab}$ 的一个估计量。此外,量 $n_{a,b}^{(t)}(i,j)$ 可以归纳地更新。事实上,

$$ n_{ab}^{t+1}(i,j) = n_{ab}^{(t)}(i,j) + \mathbb{1}\left(A_{ij}^{t} = a\right)\mathbb{1}\left(A_{ij}^{t+1} = b\right).\tag{6.13} $$

最后,初始分布也可以通过取平均来估计:

$$ \hat{\mu}^{t} = \frac{1}{\left| \left\{ (i,j) : \hat{z}_i^{(t)} = \hat{z}_j^{t} \right\} \right|} \sum_{(i,j) : \hat{z}_i^{(t)} = \hat{z}_j^{t}} A_{ij}^{t} $$

$\hat{\nu}^{t}$ 类似。这导出算法 16,用于仅知社区个数 $K$ 时对马尔可夫 SBM 聚类。注意,为节省计算时间,我们可以选择不在每个时刻更新参数。

算法 16 块相互作用参数未知时同质马尔可夫动力学的在线聚类

输入:观测图序列 $X^{1:T} = \left(X^1, \ldots, X^T\right)$;社区个数 $K$;静态图聚类算法(记为 $\mathsf{algo}$)。

输出:节点标记 $\hat{z} = \left(\hat{z}_1, \dots, \hat{z}_n\right)$。

  1. 初始化:
    • 计算 $\hat{z} \gets \mathsf{algo}\left(X^1\right)$;
    • 对 $i, j \in [N]$ 与 $a, b \in \{0,1\}$ 令 $n_{ab}(i,j) \gets 0$。

更新:

  1. for $t = 2, \cdots, T$ do
  2. 对每个节点对 $(i,j)$,用 (6.13) 更新 $n_{ab}(i,j)$;
  3. 用 (6.11) 与 (6.12) 计算 $\widehat{P}, \widehat{Q}$;
  4. 计算 $M$,使得 $M_{ij} = \sum_{a,b} n_{ab}(i,j) \log \frac{\widehat{P}_{ab}}{\widehat{Q}_{ab}}$;
  5. for $i = 1, \cdots, n$ do
  6. 对所有 $k = 1, \ldots, K$ 令 $L_{i,k} \gets \sum_{j \neq i} M_{ij}\,\mathbb{1}\left(\hat{z}_j = k\right)$;
  7. 令 $\hat{z}_i \gets \arg\max_{1 \le k \le K} L_{i,k}$。

数值结果(Numerical results)

精度随快照数的演化(Evolution of accuracy with the number of snapshots) 让我们首先用数值方法研究初始化步骤的影响。图 6.1 画出了在 50 个马尔可夫 SBM 实现上运行算法 15 所得平均精度的演化,其中初始化分别用谱聚类或随机猜测完成。显然,当谱聚类表现良好时(见图 6.1(a)),使用它比随机猜测更可取。然而引人注目的是,当初始谱聚类给出的精度很差时,似然方法能够克服这一点。例如,在图 6.1(b) 中,对第一个快照做谱聚类得到的初始聚类非常差(精度 $\approx 50\%$,因此并不比随机猜测好多少),而算法 15 确实克服了这一点,并在若干快照之后达到完美聚类。在这一特定设定下,使用谱聚类相对随机猜测没有任何优势。此外,随机猜测比谱聚类更快。

算法 15 以谱聚类或随机猜测初始化时的平均精度随快照数的演化(a:μ₁ = 2.5 log n / n,理论最小时间步数 T*theo = 14)
(a) $\mu_1 = 2.5\frac{\log n}{n}$($T_{\mathrm{theo}}^{*} = 14$)。
算法 15 以谱聚类或随机猜测初始化时的平均精度随快照数的演化(b:μ₁ = 4.0 log n / n)
(b) $\mu_1 = 4.0\frac{\log n}{n}$。
图 6.1 算法 15 在初始化分别由谱聚类或随机猜测完成时给出的平均精度的演化。合成图为具有 $n = 500$ 个节点(均分为两个簇)的马尔可夫 SBM,参数为 $\nu_1 = 1.5\frac{\log n}{n}$、$P_{11} = 0.7$ 与 $Q_{11} = 0.3$。精度对 50 次实现取平均,误差棒表示标准误。$T_{\mathrm{theo}}^{*}$ 是超过强一致阈值所需的理论最小时间步数。

未知相互作用参数(Unknown interaction parameters) 接下来,图 6.2 展示了算法 15(已知相互作用参数)与算法 16(未知相互作用参数)所得精度的对比。我们注意到,在以下所有数值实验中,我们都选择这样的稀疏设定:对单个快照做谱聚类并不比盲目随机猜测提供更多信息。虽然算法 15 提供了更好的性能(符合预期,因为它无需估计马尔可夫链转移概率),算法 16 在使用更多快照时也能达到极好的精度。

算法 15 与算法 16 在马尔可夫 SBM 上的精度对比(a:μ₁ = 0.004)
(a) $\mu_1 = 0.004$。
算法 15 与算法 16 在马尔可夫 SBM 上的精度对比(b:ν₁ = 0.016)
(b) $\nu_1 = 0.016$。
图 6.2 在线算法 15 与算法 16 在 $N = 400$、$K = 2$、$\nu_1 = 0.004$ 的马尔可夫 SBM 上所得精度的对比。结果对 25 个马尔可夫 SBM 取平均,误差棒表示标准误。

最后,我们研究算法 16 在一类马尔可夫 SBM 上的性能:其给定时间层内的期望度小于 1。这对应于一种极稀疏情形。尽管如此,如图 6.3 所示,即使 $\mu_1 = \nu_1$,只要 $P_{11} \neq Q_{11}$(见图 6.3(a)),算法 15 也表现良好。这表明,即使在最具挑战性的情形下,算法 16 也能很好地恢复社区。

极稀疏情形下算法 16 精度随快照数的演化(a:μ₁ = 0.1/N)
(a) $\mu_1 = \frac{0.1}{N}$。
极稀疏情形下算法 16 精度随快照数的演化(b:μ₁ = 0.3/N)
(b) $\mu_1 = \frac{0.3}{N}$。
图 6.3 算法 16 在极稀疏设定下所得精度随快照数的演化。我们抽取具有 $N = 300$ 个节点、两个等大小社区、参数 $\nu_1 = \frac{0.1}{N}$ 与 $P_{11} = 0.6$ 的马尔可夫 SBM。不同曲线展示 25 个马尔可夫 SBM 上的平均,误差棒对应经验标准误。

6.2.3 时序网络聚类的谱方法(Spectral Methods for Clustering Temporal Networks)

我们在 4.1 节介绍了静态图的谱方法,它们是各种组合最小化问题的松弛。其中最简单的问题是最小割(min Cut),即

$$ \operatorname{argmin}_{z \in [K]^{n}} \operatorname{Cut}(A, z), $$

其中 arg min 取遍节点集 $[n]$ 的所有可能节点标记 $z \in [K]^n$,而

$$ \operatorname{Cut}(A, z) = \sum_{i < j \,:\, z_i \neq z_j} A_{ij}. $$

现在考虑一个由其邻接矩阵列表 $\left(A^1, \ldots, A^T\right)$ 表示的时序网络。若假设各时间快照 $A^t$ 相互独立,可以简单地把经典最小割问题推广为考虑

$$ \operatorname*{argmin}_{z \in [K]^{n}} \sum_{t=1}^{T} \operatorname{Cut}\left(A^{t}, z\right). $$

由于 $\sum_{t=1}^{T} \operatorname{Cut}\left(A^{t}, z\right) = \operatorname{Cut}\left(\sum_{t=1}^{T} A^{t}, z\right)$,于是可以在时间聚合图(即以邻接矩阵 $\sum_{t=1}^{T} A^{t}$ 表示的加权图)上应用谱方法。

不幸的是,这未能考虑节点间相互作用模式中的时间相关性。举例来说,考虑这样一个网络:社区间相互作用稀疏且时间独立(因而呈尖峰状),而社区内相互作用在时间上强相关。考虑两个节点对,其相互作用模式分别为 $x_1 = (0,1,0,0,0,0,1,0,0,1)$ 与 $x_2 = (0,0,1,1,1,0,0,0,0,0)$。由于 $\|x_1\|_1 = \|x_2\|_1 = 3$,可见简单的时间聚合对时间序列 $x_1$ 与 $x_2$ 的不同时间模式不敏感,重要信息就此丢失。

一种可能的修正是把持久连接(persistent links)考虑在内。事实上,在上例中,由于 $x_1$(相应地,$x_2$)有零次(相应地,两次)$1 \to 1$ 转移,我们可以猜测 $x_2$ 来自属于同一社区的节点之间的相互作用。形式上,这可以通过考虑

$$ \operatorname*{arg\,min}_{z} \sum_{t=1}^{T} \mathrm{Cut}\left(A^{t}, z\right) + \alpha \sum_{t=2}^{T} \mathrm{PerCut}\left(A^{t-1}, A^{t}, z\right), $$

来实现,其中

$$ \mathrm{PerCut}\left(A^{t-1}, A^{t}, z\right) = \sum_{i, j \,:\, z_i \neq z_j} A_{ij}^{t-1} A_{ij}^{t} $$

计算从时刻 $t-1$ 到时刻 $t$ 割中持久连接的数目。我们进一步注意到

$$ \mathrm{PerCut}\left(A^{t-1}, A^{t}, z\right) = \mathrm{Cut}\left(A^{t-1} \odot A^{t}, z\right) $$

其中 $\odot$ 表示矩阵的逐元乘积。下一小节将为在时序网络聚类中考虑持久边的直觉提供论证。

带马尔可夫边动力学的度校正时序 SBM(Degree-corrected temporal SBM with Markov edge dynamics) 我们首先介绍马尔可夫 SBM 的一个度校正版本。一个具有 $n$ 个节点、$K$ 个块和 $T$ 个快照的度校正时序随机分块模型可以由对称邻接张量 $A \in \{0,1\}^{n\times n\times T}$(对角元为零)的概率分布

$$ \mathbb{P}(A \mid Z, F, \theta) = \prod_{1 \leq i < j \leq n} F_{z_i z_j}^{\theta_i \theta_j}\left(A_{ij}^{1}, \ldots, A_{ij}^{T}\right)\tag{6.14} $$

描述,其中 $z = (z_1, \ldots, z_n)$ 是社区指派($z_i \in [K]$ 指示节点 $i$ 的社区),$F = \left(F_{k\ell}^{xy}\right)$ 是 $\{0,1\}^T$ 上的一族概率分布,$\theta = (\theta_1, \ldots, \theta_n)$ 是节点特定的度校正参数向量,满足 $0 \leq \theta_i < \infty$。

在下文中,我们限于具有马尔可夫边动力学的同质块间相互作用情形:节点的静态社区标记从所有节点标记的集合 $[K]$ 中均匀随机抽取,且

$$ F_{z_i z_j}^{\theta_i \theta_j}(x) = \left\{ \begin{array}{ll} \mu_{x_1}^{\theta_i \theta_j} \prod_{t=2}^{T} P_{x_{t-1}, x_t}^{\theta_i \theta_j} & \quad \text{若 } z_i = z_j, \\ \nu_{x_1}^{\theta_i \theta_j} \prod_{t=2}^{T} Q_{x_{t-1}, x_t}^{\theta_i \theta_j} & \quad \text{其他}, \end{array} \right.\tag{6.15} $$

其中初始分布为

$$ \mu^{\theta_i \theta_j} = \binom{1 - \theta_i \theta_j \mu_1}{\theta_i \theta_j \mu_1}, \quad \nu^{\theta_i \theta_j} = \binom{1 - \theta_i \theta_j \nu_1}{\theta_i \theta_j \nu_1}, $$

转移概率矩阵为

$$ 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}, \quad 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}. $$

参数 $\theta_i$($i = 1, \ldots, n$)刻画了某些节点比其他节点更倾向于发起新连接的事实,与度校正分块模型(Karrer and Newman, 2011)类似。为保持模型简单,我们不在 $P_{11}$ 前加度校正参数;因此一旦连接建立,保持其活跃的概率就是 $P_{11}$ 或 $Q_{11}$。此外,我们假设 $\min_{i,j}\{\theta_i \theta_j \delta\} \le 1$,其中 $\delta = \max\{\mu_1, \nu_1, P_{01}, Q_{01}\}$。最后,我们对度校正参数做归一化,使得对所有 $k$ 有 $\sum_i \mathbb{1}(z_i = k)\theta_i = \sum_i \mathbb{1}(z_i = k)$。

最大似然估计量(Maximum likelihood estimator)

命题 6.2 度校正马尔可夫 SBM 的最大似然估计量(Avrachenkov et al., 2021b)

令 $\rho_a^{\theta_i \theta_j} = \log \dfrac{\mu_a^{\theta_i \theta_j}}{\nu_a^{\theta_i \theta_j}}$,$\ell_{ab}^{\theta_i \theta_j} = \log \dfrac{P_{ab}^{\theta_i \theta_j}}{Q_{ab}^{\theta_i \theta_j}} - \log \dfrac{P_{00}^{\theta_i \theta_j}}{Q_{00}^{\theta_i \theta_j}}$。由 (6.14)–(6.15) 定义的度校正马尔可夫 SBM 的一个最大似然估计量,是任意使下式最大化的社区指派 $\hat{z}$:

$$ \sum_{\substack{i,j \\ z_i = z_j}} \left\{ A_{ij}^{1}\left(\rho_1^{\theta_i \theta_j} - \rho_0^{\theta_i \theta_j}\right) + \rho_0^{\theta_i \theta_j} + \left(A_{ij}^{1} - A_{ij}^{T}\right)\ell_{10}^{\theta_i \theta_j} \right\} + \sum_{\substack{i,j \\ z_i = z_j}} \sum_{t=2}^{T} \left\{ \left(\ell_{01}^{\theta_i \theta_j} + \ell_{10}^{\theta_i \theta_j}\right)\left(A_{ij}^{t} - A_{ij}^{t-1}A_{ij}^{t}\right) + \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\} $$

(最大化取遍所有社区指派 $z \in [K]^n$。)

证明 命题 6.2

由时间马尔可夫性,模型的对数似然可以写成 $\log \mathbb{P}(A \mid z, \theta) = \log \mathbb{P}(A^1 \mid z, \theta) + \sum_{t=2}^{T} \mathbb{P}(A^t \mid A^{t-1}, z, \theta)$。(校勘:右端第二项原书漏印 $\log$,应为 $\sum_{t=2}^T \log \mathbb{P}(A^t \mid A^{t-1}, z, \theta)$。)记 $\rho_a^{\theta_i \theta_j} = \log \frac{\mu_a^{\theta_i \theta_j}}{\nu_a^{\theta_i \theta_j}}$,我们得到

$$ \begin{array}{l} \displaystyle \log \mathbb{P}(A^1 \mid z, \theta) = \frac{1}{2}\sum_{i,j}\sum_{a} \delta(A_{ij}^{1}, a)\Big( \delta(z_i, z_j)\rho_a^{\theta_i \theta_j} + \log \nu_a^{\theta_i \theta_j} \Big) \\ \displaystyle \quad = \frac{1}{2}\sum_{i,j} \delta(z_i, z_j)\sum_{a} \delta(A_{ij}^{1}, a)\rho_a^{\theta_i \theta_j} + c_1(A), \end{array} $$

其中 $c_1(A) = \frac{1}{2}\sum_{i,j}\sum_{a} \delta(A_{ij}^{1}, a)\log \nu_a^{\theta_i \theta_j}$ 不依赖于社区结构。类似地,记 $R_{ab}^{\theta_i \theta_j} = \log \frac{P_{ab}^{\theta_i \theta_j}}{Q_{ab}^{\theta_i \theta_j}}$,我们得到 $\log \mathbb{P}(A^t \mid A^{t-1}, z, \theta)$ 等于

$$ \begin{array}{rl} \displaystyle \frac{1}{2}\sum_{i,j}\sum_{a,b} \delta(A_{ij}^{t-1}, a)\delta(A_{ij}^{t}, b)\Big( \delta(z_i, z_j)R_{ab}^{\theta_i \theta_j} + \log Q_{ab}^{\theta_i \theta_j} \Big) & \\ = \displaystyle \frac{1}{2}\sum_{i,j} \delta(z_i, z_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{array} $$

其中 $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}$ 不依赖于社区结构。简单的计算表明

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

且 $\sum_{a,b} \delta(A_{ij}^{t-1}, a)\delta(A_{ij}^{t}, b)R_{ab}^{\theta_i \theta_j}$ 等于

$$ \begin{array}{l} R_{00}^{\theta_i \theta_j} + A_{ij}^{t-1}\left(R_{10}^{\theta_i \theta_j} - R_{00}^{\theta_i \theta_j}\right) + A_{ij}^{t}\left(R_{01}^{\theta_i \theta_j} - R_{00}^{\theta_i \theta_j}\right) \\ \qquad\qquad + A_{ij}^{t-1}A_{ij}^{t}\left(R_{11}^{\theta_i \theta_j} - R_{01}^{\theta_i \theta_j} - R_{10}^{\theta_i \theta_j} + R_{00}^{\theta_i \theta_j}\right) \\ = 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} + A_{ij}^{t-1}A_{ij}^{t}\left(\ell_{11}^{\theta_i \theta_j} - \ell_{01}^{\theta_i \theta_j} - \ell_{10}^{\theta_i \theta_j}\right). \end{array} $$ 查看学习笔记对这一“简单的计算”的逐步展开

汇总上述观察,我们得到 $\log \mathbb{P}(A \mid z, \theta)$ 等于

$$ \begin{array}{l} \displaystyle c(A) + \frac{1}{2}\sum_{\substack{i,j \\ z_i = z_j}} \left\{ A_{ij}^{1}\left(\rho_1^{\theta_i \theta_j} - \rho_0^{\theta_i \theta_j}\right) + \rho_0^{\theta_i \theta_j} + \left(A_{ij}^{1} - A_{ij}^{T}\right)\ell_{10}^{\theta_i \theta_j} \right\} \\ \displaystyle \quad + \frac{1}{2}\sum_{\substack{i,j \\ z_i = z_j}} \sum_{t=2}^{T} \left\{ \left(\ell_{01}^{\theta_i \theta_j} + \ell_{10}^{\theta_i \theta_j}\right)\left(A_{ij}^{t} - A_{ij}^{t-1}A_{ij}^{t}\right) + \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\}, \end{array} $$

其中 $c(A) = \sum_t c_t(A)$ 不依赖于 $z$。命题由此得证。

查看学习笔记完整证明(含跳步整合与校勘)

命题 6.2 推导的 MLE 比独立地把所有快照求和得到的 MLE 更复杂。特别地,项 $A_{ij}^{t-1}A_{ij}^{t}$ 刻画了两个相邻快照之间的持久边。记 $A_{\mathrm{pers}}^{t} = A^{t-1} \odot A^{t}$ 为邻接矩阵 $A^{t-1}$ 与 $A^{t}$ 的逐元乘积。则 $A_{\mathrm{pers}}^{t}$ 是包含 $t-1$ 与 $t$ 之间持久边的图的邻接矩阵,而 $A_{\mathrm{new}}^{t} = A^{t} - A_{\mathrm{pers}}^{t}$ 对应于包含时刻 $t-1$ 到时刻 $t$ 之间新生出现的边的图。

假设快照数 $T$ 很大,我们可以忽略边界项,命题 6.2 表达的 MLE 就化为最大化

$$ \sum_{t=2}^{T} \sum_{\substack{i,j \\ z_i = z_j}} \left( \left(\ell_{01}^{\theta_i \theta_j} + \ell_{10}^{\theta_i \theta_j}\right)\left(A_{ij}^{t} - A_{ij}^{t-1}A_{ij}^{t}\right) + \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). $$

这个表达式可以进一步化简,表示为正则化模块度。回忆:给定加权图 $W$、划分 $z$ 与分辨率参数 $\gamma$,正则化模块度定义为(见 4.2 节与式 (4.21))

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

其中 $d_i = \sum_j W_{ij}$,$m = \frac{1}{2}\sum_i d_i$。

引理 6.1 稀疏设定下 MLE 近似为正则化模块度最大化

假设 $P^{\theta_i \theta_j}$ 与 $Q^{\theta_i \theta_j}$ 非退化,且 $\mu^{\theta_i \theta_j}$(相应地,$\nu^{\theta_i \theta_j}$)是 $P^{\theta_i \theta_j}$(相应地,$Q^{\theta_i \theta_j}$)的平稳分布。在 $P_{01}$ 与 $Q_{01}$ 很小的稀疏设定下,MLE 近似地最大化 $\mathcal{M}(W, z, \gamma)$,其中 $W$ 定义为

$$ W = \sum_{t=2}^{T} \left( \alpha A_{\mathrm{new}}^{t} + \beta A_{\mathrm{pers}}^{t} \right),\tag{6.16} $$

其中

$$ \alpha = \log \frac{P_{01}}{Q_{01}} + \log \frac{1 - P_{11}}{1 - Q_{11}}, \quad \text{且} \quad \beta = \log \frac{P_{11}}{Q_{11}},\tag{6.17} $$ $$ \text{且} \quad \gamma = (P_{01} - Q_{01})\,\frac{\alpha\left(\mu_1 + (K-1)\nu_1\right) + (\beta - \alpha)\left(\mu_1 P_{11} + (K-1)\nu_1 Q_{11}\right)}{K}. $$
Tips:这是 6.2.3 节的桥梁结果:马尔可夫边动力学下的 MLE 近似化为“新生边加权 $\alpha$、持久边加权 $\beta$”的加权图上的正则化模块度最大化,从而可用第 4 章的归一化谱聚类求解(算法 17)。原书给出的两个 $\gamma$ 版本都不闭合;对式 (6.18) 先完成时间求和后,校勘值为 $\gamma=K(P_{01}-Q_{01})/S$,见证明后的提示。
证明 引理 6.1

由于 $P_{01}, Q_{01} = o(1)$,一阶 Taylor 展开给出

$$ \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\left(P_{01}^{2} + Q_{01}^{2}\right), $$

以及 $\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}}$、$\ell_{11}^{\theta_i \theta_j} \approx \log \frac{P_{11}}{Q_{11}}$。把这些近似代入 MLE 表达式,得到最大化

$$ \sum_{t=2}^{T} \sum_{i,j} \delta(z_i, z_j)\left( \tilde{a}_{ij}^{t} - \theta_i \theta_j\left(P_{01} - Q_{01}\right) \right)\tag{6.18} $$

其中 $\tilde{a}_{ij}^{t} = \alpha\left(A_{\mathrm{new}}^{t}\right)_{ij} + \beta\left(A_{\mathrm{pers}}^{t}\right)_{ij}$。由于 $\mu$ 与 $\nu$ 是平稳分布,

$$ \mathbb{E}\left(A_{\mathrm{new}}^{t}\right)_{ij} = \left\{ \begin{array}{ll} \theta_i \theta_j \mu_1(1 - P_{11}) & \text{若 } z_i = z_j, \\ \theta_i \theta_j \nu_1(1 - Q_{11}) & \text{其他}, \end{array} \right. $$ $$ \mathbb{E}\left(A_{\mathrm{pers}}^{t}\right)_{ij} = \left\{ \begin{array}{ll} \theta_i \theta_j \mu_1 P_{11} & \text{若 } z_i = z_j, \\ \theta_i \theta_j \nu_1 Q_{11} & \text{其他}. \end{array} \right. $$

因此,利用 $W_{ij} = \sum_{t=2}^{T} \tilde{a}_{ij}^{t}$,我们得到

$$ \mathbb{E}W_{ij} = \left\{ \begin{array}{ll} (T-1)\theta_i \theta_j \mu_1\left(\alpha(1 - P_{11}) + \beta P_{11}\right) & \text{若 } z_i = z_j, \\ (T-1)\theta_i \theta_j \nu_1\left(\alpha(1 - Q_{11}) + \beta Q_{11}\right) & \text{其他}. \end{array} \right. $$

由于社区标记均匀随机抽取,且 $\theta_i$ 已适当归一化,期望度 $\bar{d}_i$ 等于

$$ (T-1)\theta_i n\,\frac{\mu_1\left(\alpha(1 - P_{11}) + \beta P_{11}\right) + (K-1)\nu_1\left(\alpha(1 - Q_{11}) + \beta Q_{11}\right)}{K}, $$

并有 $\bar{m} = \dfrac{n^2}{2}\,\dfrac{\mu_1\left(\alpha(1 - P_{11}) + \beta P_{11}\right) + (K-1)\nu_1\left(\alpha(1 - Q_{11}) + \beta Q_{11}\right)}{K}$。因此,我们观察到 $\theta_i \theta_j(P_{01} - Q_{01}) = \gamma\dfrac{\bar{d}_i \bar{d}_j}{2\bar{m}}$,其中 $\gamma = (P_{01} - Q_{01})(T-1)\,\dfrac{\mu_1\left(\alpha(1 - P_{11}) + \beta P_{11}\right) + (K-1)\nu_1\left(\alpha(1 - Q_{11}) + \beta Q_{11}\right)}{K}$。利用式 (6.18) 结束证明。

查看学习笔记完整证明(含跳步整合与校勘)

结合新生边与持久边的时序谱聚类(Temporal spectral clustering combining new and persistent edges) 根据上一小节的分析,MLE 近似地由下式的解给出:

$$ \operatorname*{arg\,max}_{z \in [K]^{n}} \mathcal{M}(W, z, \gamma), $$

其中 $W$ 由式 (6.16) 定义,$\gamma$ 是适当的分辨率参数。这个优化问题一般是 NP 完全的(Brandes et al., 2007),但可以用连续松弛近似求解。我们可以选取松弛,使得该优化问题化为加权图 $W$ 上的归一化谱聚类算法(见 4.4.2 节)。我们注意到,为计算 $W$ 的归一化拉普拉斯矩阵,应当限制 $\alpha, \beta \geq 0$,而公式 (6.17) 并不必然保证这一点。我们将此总结为算法 17。

算法 17 马尔可夫边动力学与静态节点标记下时序网络的谱聚类

输入:邻接矩阵 $A^1, \cdots, A^T$;簇数 $K$;参数 $\alpha, \beta$。

输出:预测的社区标签 $\hat{z} \in [K]^N$。

过程:

  • 令 $W = \sum_{t=2}^{T}\left(\alpha A_{\mathrm{new}}^{t} + \beta A_{\mathrm{pers}}^{t}\right)$,其中 $A_{\mathrm{new}}^{t} = A^{t} - A^{t-1} \odot A^{t}$、$A_{\mathrm{pers}}^{t} = A^{t-1} \odot A^{t}$;
  • 计算 $\mathcal{L} = I_n - D^{-1/2}WD^{-1/2}$,其中 $D = \mathrm{diag}(W\mathbb{1}_n)$;
  • 计算 $\widehat{X} \in \mathbb{R}^{N \times K}$,其各列为 $\mathcal{L}$ 对应于 $K$ 个最小特征值的 $K$ 个标准正交特征向量。

返回:$\hat{z} \gets k\text{-means}\left(D^{-1/2}\widehat{X}, K\right)$。

数值结果(Numerical results)

合成数据(Synthetic data) 我们首先考察算法 17 中参数 $\alpha$ 与 $\beta$ 选择的影响。为此,令 $\alpha = 1$,并在图 6.4 中对不同的 $\beta$ 画出在 25 个马尔可夫边动力学随机分块模型实现上得到的平均精度。虽然在时间聚合图上做谱聚类(对应于 $\beta = 1$)效果不错,但引人注目的是,其他 $\beta$ 值给出了甚至更好的结果。$\beta$ 的选择取决于持久相互作用的概率。例如,若 $P_{11} > Q_{11}$(图 6.4(a)),则宜取 $\beta > 1$;而若 $P_{11} < Q_{11}$(图 6.4(b)),取大的 $\beta$ 会受到惩罚。这与公式 (6.17) 推荐的 $\alpha$、$\beta$ 取值一致。

算法 17 在合成时序 SBM 上不同 β 取值下的精度(a:P₁₁ = 0.9,同社区边持久性更高)
(a) $P_{11} = 0.9$。
算法 17 在合成时序 SBM 上不同 β 取值下的精度(b:P₁₁ = 0.1,同社区边持久性低于跨社区)
(b) $P_{11} = 0.1$。
图 6.4 算法 17 在具有 300 个节点、$K = 3$ 个块、平稳马尔可夫边演化 $\mu_1 = 0.04$、$\nu_1 = 0.02$ 与 $Q_{11} = 0.3$ 的时序 SBM 上的精度。结果对 25 个合成图实现取平均,误差棒表示标准差。

高中学生的社交网络(Social networks of high school students) 我们考察连续三年从法国马赛 Lycée Thiers 高中采集的三个数据集(Fournet and Barrat, 2014; Mastrandrea et al., 2015)。我们在引言中介绍过这些数据集。特别地,节点对应学生,相互作用对应近距离接触事件,社区对应班级,各维度见表 1.1。

我们假设各年相互作用的时间特征相似。然后用 2011 年数据集估计转移概率矩阵 $P$ 与 $Q$,并把它们用于 2012 与 2013 年数据集的聚类。我们假设 $\theta_i = 1$(无度校正)。马尔可夫链转移概率矩阵的标准估计量(Billingsley, 1961)给出

$$ \widehat{P} = \begin{pmatrix} 0.9992 & 0.0008 \\ 0.37 & 0.63 \end{pmatrix} \quad \text{且} \quad \widehat{Q} = \begin{pmatrix} 0.999967 & 3.3 \times 10^{-5} \\ 0.48 & 0.52 \end{pmatrix}. $$

利用 (6.17),得到 $\hat{\alpha} = 2.9$ 与 $\hat{\beta} = 0.18$。我们在图 6.5(b) 中观察到,在 2013 年数据集上,这组参数比简单地在时间聚合图上做谱聚类($\alpha = \beta = 1$)给出更好的精度。对 2012 年数据集(图 6.5(a)),这一改进不那么明显可见。

算法 17 在 2012 年高中数据集上的精度:均匀权重 (1,1)(蓝)与用 2011 年数据估计的调整权重 (2.9, 0.18)(橙)对比
(a) 2012 年。
算法 17 在 2013 年高中数据集上的精度:调整权重(橙)明显优于均匀权重(蓝)
(b) 2013 年。
图 6.5 算法 17 在 2012 与 2013 年高中数据集上的精度,分别使用均匀权重 $\alpha = \beta = 1$(蓝色)与调整后的 $\alpha, \beta$(其取值用 2011 年数据预测,橙色)。

为理解算法 17 为什么在 2013 年比 2012 年表现更好,我们在表 6.1 中列出了对每个数据集单独估计的时间转移概率与聚类权重 $\hat{\alpha}, \hat{\beta}$。对 2012 年,社区内边持久性 $\widehat{P}_{11}$ 与社区间边持久性 $\widehat{Q}_{11}$ 之差很小,这意味着持久边没有为区分社区带来太多额外信息($\hat{\beta} \approx 0$)。对 2011 与 2013 年,这一差异更大,表明边持久性含有可用于以更高精度恢复社区的信息。

数据集 $\widehat{P}_{01}$ $\widehat{Q}_{01}$ $\widehat{P}_{11}$ $\widehat{Q}_{11}$ $\widehat{\alpha}$ $\widehat{\beta}$ $\widehat{\beta}/\widehat{\alpha}$
20110.000800.0000330.630.522.90.580.060
20120.000500.0000110.570.563.80.010.003
20130.001500.0000140.640.404.50.070.015
表 6.1 对每个数据集单独估计的马尔可夫链转移概率与调整后的聚类权重。

6.2.4 用经验转移概率进行长时间跨度聚类(原文章节标题作 Clustering for Long Time Horizon Using Empirical Transition Rates)

我们继续研究例 6.2 定义的具有静态成员与同质马尔可夫相互作用核的时序 SBM。记 $P, Q$ 为转移概率矩阵。让我们考虑快照数 $T$ 趋于无穷而 $N$ 保持有界的情形。主要思想是利用马尔可夫链的遍历性,用标准技术估计参数,然后做推断。目前我们假设相互作用参数 $P, Q$ 已知,但 $K$ 未知。$P, Q$ 同样未知的情形见注 6.3。

回忆公式 (6.10) 给出了 $P(i,j)$——节点对 $(i,j)$ 的模式相互作用演化的转移概率矩阵——的一致估计量。于是,一旦所有 $P(i,j)$ 都以良好精度已知,我们就可以利用对 $P, Q$ 的了解来判别节点 $i$ 与 $j$ 是否在同一块中,并用这些数据在节点集上构造一个相似图。这导出算法 18:它不需要先验地知道块数,而是把块数作为副产品估计出来。注意,该算法是为同质相互作用阵列量身定制的。

算法 18 基于经验转移概率的聚类

输入:观测相互作用张量 $\left(A_{ij}^{t}\right)$;转移概率矩阵 $P, Q$。

输出:估计的节点标记 $\hat{z} = \left(\hat{z}_1, \dots, \hat{z}_n\right)$;估计的社区个数 $\widehat{K}$。

  1. 令 $V \gets \{1, \ldots, n\}$,$E \gets \varnothing$。
  2. for 所有无序节点对 $ij$ do
  3. 用 (6.10) 对 $a, b = 0, 1$ 计算 $\widehat{P}_{ab}(i,j)$。
  4. if 对某个 $a, b$ 有 $\left| \widehat{P}_{ab}(i,j) - P_{ab} \right| \le \frac{1}{2}\left| P_{ab} - Q_{ab} \right|$ then
  5. 令 $E \gets E \cup \{ij\}$。
  6. 计算 $\mathcal{C} \gets$ 图 $G = (V, E)$ 中的连通分量集合,令 $\widehat{K} \gets |\mathcal{C}|$,并令 $(C_1, \dots, C_{\widehat{K}}) \gets$ 按任意顺序列出的 $\mathcal{C}$ 的成员。
  7. for $i = 1, \ldots, n$ do
  8. 令 $\hat{z}_i \gets$ 满足 $C_k \ni i$ 的唯一 $k$。
定理 6.2 长时间跨度下算法 18 的完全恢复

考虑一个具有 $n$ 个节点、$K$ 个社区和 $T$ 个快照的同质马尔可夫 SBM。假设 $n$ 固定,且转移概率矩阵 $P, Q$ 已知。则当 $T$ 趋于无穷时,只要演化不是静态的且 $P \neq Q$,算法 18 以高概率把每个节点正确分类。

Tips:这是 $n$ 固定、$T \to \infty$ 这一与第 4 章互补的渐近情形:遍历性使每个节点对的经验转移概率一致收敛,把“是否同社区”的判别化为相似图上的连通分量计算,且无需预知社区数 $K$。
证明 定理 6.2

对 $a, b \in \{0,1\}$,令 $n_a(i,j) = \sum_b n_{ab}(i,j)$,其中 $n_{ab}(i,j)$ 记录节点对 $(i,j)$ 之间观测到的 $a \to b$ 转移次数。随机变量

$$ \xi_{ab}(i,j) = \frac{n_{ab}(i,j) - n_a(i,j)P_{ab}(i,j)}{\sqrt{n_a(i,j)}} $$

的分布趋于一个零均值、有限方差的正态分布,其方差由 $\lambda_{(ab),(cd)} = \delta_{ac}\left(\delta_{bd}P_{ab}(i,j) - P_{ab}(i,j)P_{a,d}(i,j)\right)$ 给出(见 Billingsley, 1961, Theorem 3.1 与公式 (3.13))。因此,对任意 $\alpha > 0$,

$$ \mathbb{P}\Big( \left| \widehat{P}_{ab}(i,j) - P_{ab}(i,j) \right| \geq \alpha \Big) = \mathbb{P}\Big( \left| \xi_{ab}(i,j) \right| \geq \alpha \sqrt{n_a(i,j)} \Big),\tag{6.19} $$

且该量随 $T$ 趋于无穷而趋于零。

由模型的可辨识性,$P \neq Q$。因此不失一般性,我们可以假设 $P_{01} \neq Q_{01}$,并选取 $\alpha$ 使得 $0 < \alpha < \frac{P_{01} - Q_{01}}{2}$。当 $\widehat{P}_{01}(i,j) > \frac{P_{01} + Q_{01}}{2}$ 时,节点 $i$ 与 $j$ 被预测为在同一社区,犯错的概率为

$$ \mathbb{P}\left( \left| \widehat{P}_{01}(i,j) - P_{01}(i,j) \right| \geq \alpha \right). $$

由并集界,所有节点都被正确分类的概率以

$$ \frac{n(n-1)}{2}\max_{ij} \mathbb{P}\left( \left| \widehat{P}_{01}(i,j) - P_{01}(i,j) \right| \geq \alpha \right), $$

为界,其中最大值取遍所有节点对 $ij$。由式 (6.19),对所有节点对 $ij$,有 $\mathbb{P}\Big( \left| \widehat{P}_{01}(i,j) - P_{01}(i,j) \right| \geq \alpha \Big) \to 0$。因此,当 $T \to \infty$ 时,所有节点都几乎必然被正确分类。

查看学习笔记的条件补全与印刷算法审计
注 6.3 $P, Q$ 未知时的处理

若 $P$ 与 $Q$ 未知,我们可以增加一步:把估计的转移矩阵 $\widehat{P}(i,j)$ 聚成两类(例如使用 k-means)。

6.3 社区成员的马尔可夫演化(Markovian Evolution of Community Memberships)

本节关注这样的时序网络聚类:其成员结构服从马尔可夫链,但相互作用结构是时间独立的。具体地,记 $z_{it} \in [K]$ 为节点 $i$ 在时刻 $t$ 的组成员。那么,跨节点来看,随机变量 $\left(z_{it}\right)_{1 \le t \le T}$ 独立同分布。对每个节点 $i$,组成员 $z_{i\cdot} = (z_{i1}, \cdots, z_{iT})$ 服从一个不可约且非周期的马尔可夫链,由

$$ \mathbb{P}\left(z_{i\cdot}\right) = \alpha_{z_{i1}} \prod_{t=2}^{T} \pi_{z_{i,t-1}, z_{it}},\tag{6.20} $$

给出,其中 $\alpha$ 是初始分布,$\pi$ 是转移概率矩阵。在节点标签的条件下,各边相互独立,且对所有 $i < j$ 与所有 $t$,有

$$ A_{ij}^{t} \mid z_{it}, z_{jt} \sim \mathrm{Ber}\left(p_{z_{it} z_{jt}}\right). $$

因此,邻接矩阵序列 $A^{1:T} = \left(A^1, \cdots, A^T\right)$ 的似然为

$$ \mathbb{P}\left(A^{1:T} \mid Z\right) = \prod_{i=1}^{n} \alpha\left(z_{i1}\right) \prod_{t=2}^{T} \pi_{z_{i,t-1}, z_{it}} \prod_{t=1}^{T} \prod_{i < j} p_{z_{it} z_{jt}}^{A_{ij}^{t}}\left(1 - p_{z_{it} z_{jt}}\right)^{1 - A_{ij}^{t}}. $$

6.3.1 变分期望最大化算法(Variational Expectation–Maximization Algorithm)

让我们首先假设 $K$ 已知,且 $\alpha$ 是 $\pi$ 的平稳分布。我们的目标是估计组成员 $Z = \left(z_{it}\right)_{1 \le i \le n, 1 \le t \le T}$ 以及模型参数 $\theta = (\pi, P)$,其中 $P = \left(p_{k\ell}\right)_{k,\ell\in[K]}$。

当 $n$ 或 $T$ 很大时,似然的全局最大化是不可行的,而期望最大化(EM)算法(Dempster et al., 1977)提供了一种寻找局部最大值的方法。EM 算法计算给定观测 $A^{1:T}$ 时 $Z$ 的条件分布。然而,在我们的情形中,由于依赖性,该分布不能分解为 $n$ 个节点上的乘积。事实上,我们有

$$ \mathbb{P}\left(Z \mid A^{1:T}\right) = \mathbb{P}\left(z_{\cdot 1} \mid A^{1}\right) \prod_{t=2}^{T} \mathbb{P}\left(z_{\cdot t} \mid z_{\cdot t-1}, A^{t}\right), $$

其中 $z_{\cdot t} = (z_{1t}, \cdots, z_{nt})$ 表示时刻 $t$ 的社区标签。不幸的是,分布 $\mathbb{P}\left(Z^{t} \mid Z^{t-1}, A^{t}\right)$ 无法进一步分解,因为随机变量 $z_{it} \mid A_{ij}^{t}$ 与 $z_{jt} \mid A_{ij}^{t}$ 并不独立。事实上,在时刻 $t$ 观测到 $i$ 与 $j$ 之间的一条边会增大 $z_{it} = z_{jt}$ 的可能性。变分近似引入如下一类概率分布 $\mathbb{Q}$:

$$ \mathbb{Q}_{\tau}(Z) = \prod_{i=1}^{n} \mathbb{Q}_{\tau}\left(z_{i\cdot}\right) = \prod_{i=1}^{n} \mathbb{Q}_{\tau}\left(z_{i1}\right) \prod_{t=2}^{T} \mathbb{Q}_{\tau}\left(z_{it} \mid z_{it-1}\right). $$

我们引入 $\tau(i,k) = \mathbb{Q}\left(z_{i1} = k\right)$ 与 $\tau(t,i,k,\ell) = \mathbb{Q}\left(z_{it} = \ell \mid z_{it-1} = k\right)$。于是,在 $\mathbb{Q}$ 下,$(z_{i1}, \dots, z_{iT})$ 的分布是一个时间非齐次的马尔可夫链,其转移为 $\tau(t,i,k,\ell)$、初始分布为 $\tau(i,k)$。特别地,$\sum_{k=1}^{K} \tau(i,k) = 1$、$\sum_{\ell=1}^{K} \tau(t,i,k,\ell) = 1$,且

$$ \mathbb{Q}(Z) = \prod_{i=1}^{n} \prod_{k=1}^{K} \tau(i,k)^{\mathbb{1}(z_{it} = k)} \prod_{t=2}^{T} \prod_{1 \le k, \ell, K} \tau(t,i,k,\ell)^{\mathbb{1}(z_{it-1} = k)\mathbb{1}(z_{it} = \ell)}. $$

边际分布 $\tau_{\mathrm{marg}}(t,i,k) = \mathbb{Q}\left(z_{it} = k\right)$ 由

$$ \begin{array}{l} \displaystyle \tau_{\mathrm{marg}}(1, i, k) = \tau(i,k), \\ \displaystyle \tau_{\mathrm{marg}}(t, i, k) = \sum_{\ell=1}^{K} \tau_{\mathrm{marg}}(t-1, i, \ell)\,\tau(t, i, \ell, k) \end{array} $$

递归计算。

变分期望最大化(VEM)算法(Matias and Miele, 2017)则寻求最大化

$$ J(\theta, \tau) = \mathbb{E}_{\mathbb{Q}}\left( \log \mathbb{P}\left(A^{1:T}, Z\right) \right) + \mathcal{H}(\mathbb{Q}), $$

其中 $\mathcal{H}(\mathbb{Q})$ 表示 $\mathbb{Q}$ 的熵。因此,$J(\theta, \tau)$ 等于

$$ \begin{array}{l} \displaystyle \sum_{i=1}^{n} \sum_{k=1}^{K} \tau(i,k)\big[ \log \alpha_k - \log \tau(i,k) \big] \\ \displaystyle \quad + \sum_{t=2}^{T} \sum_{i=1}^{n} \sum_{1 \le k, \ell \le K} \tau_{\mathrm{marg}}(t-1, i, k)\,\tau(t, i, k, \ell) \times \big[ \log \pi_{k\ell} - \log \tau(t,i,k,\ell) \big] \\ \displaystyle \quad + \sum_{t=1}^{T} \sum_{1 \le i < j \le n} \sum_{1 \le k, \ell \le K} \tau_{\mathrm{marg}}(t, i, k)\,\tau_{\mathrm{marg}}(t, j, \ell) \times \log \left( \mathrm{Ber}\left(p_{z_{it} z_{jt}}\right)\left(A_{ij}^{t}\right) \right), \end{array} $$

其中

$$ \mathrm{Ber}\left(p_{z_{it} z_{jt}}\right)\left(A_{ij}^{t}\right) = \left\{ \begin{array}{ll} p_{z_{it} z_{jt}} & \text{若 } A_{ij}^{t} = 1, \\ 1 - p_{z_{it} z_{jt}} & \text{其他}. \end{array} \right. $$

优化迭代地进行。在第 $k$ 步,给定当前估计 $\left(\tau^{k}, \theta^{k}\right)$,我们执行下面两个子步骤:

  1. VE 步: 计算 $\tau^{k+1} = \arg\max_{\tau} J\left(\theta^{k}, \tau\right)$;
  2. M 步: 计算 $\theta^{k+1} = \arg\max_{\theta} J\left(\theta, \tau^{k+1}\right)$。

下述引理给出更新 $\tau^{k+1}$ 与 $\theta^{k+1}$ 的值。

引理 6.3 VEM 更新的显式形式

值 $\widehat{\tau} = \arg\max_{\tau} J(\theta, \tau)$ 满足

$$ \widehat{\tau}(t,i,k,\ell) \propto \pi_{k\ell} \prod_{j=1}^{n} \prod_{k'=1}^{K} \left( \mathrm{Ber}\left(p_{z_{it} z_{jt}}\right)\left(A_{ij}^{t}\right) \right)^{\widehat{\tau}_{\mathrm{marg}}(t,j,k')}, $$

其中比例关系保证 $\tau$ 的归一化约束。类似地,$\widehat{\theta} = \arg\max_{\theta} J(\tau, \theta)$ 由 $\widehat{\theta} = \left(\widehat{\pi}, \widehat{P}\right)$ 给出,使得

$$ \begin{array}{rl} & \displaystyle \widehat{\pi}_{k\ell} \propto \sum_{t=2}^{T} \sum_{i=1}^{n} \tau_{\mathrm{marg}}(t-1, i, k)\,\tau(t,i,k,\ell), \\ & \displaystyle \widehat{p}_{k\ell} = \frac{\sum_{t=1}^{T} \sum_{1 \le i,j \le n} \tau_{\mathrm{marg}}(t,i,k)\,\tau_{\mathrm{marg}}(t,i,\ell)\,\mathbb{1}\left(A_{ij}^{t} \neq 0\right)}{\sum_{t=1}^{T} \sum_{1 \le i,j \le n} \tau_{\mathrm{marg}}(t,i,k)\,\tau_{\mathrm{marg}}(t,i,\ell)}. \end{array} $$
证明 引理 6.3

证明由对 $J(\tau, \theta)$ 直接求导得到。例如,我们有

$$ \begin{array}{rl} & \displaystyle \frac{\partial J}{\partial \tau(t,i,k,\ell)} = \tau_{\mathrm{marg}}(t-1,i,k)\big[ \log \pi_{k\ell} - \log \tau(t,i,k,\ell) + 1 \big] \\ & \displaystyle \qquad + \tau(t,i,k,\ell)\,\tau_{\mathrm{marg}}(t-1,i,k) \log \left( \mathrm{Ber}\left(p_{z_{it} z_{jt}}\right)\left(A_{ij}^{t}\right) \right) \end{array} $$

令该导数为零,即得到所述的 $\widehat{\tau}(t,i,k,\ell)$ 表达式。(校勘:上式第二项中 $\mathrm{Ber}$ 的下标 OCR 误植为 $z_i z_j t$,原书作 $p_{z_{it}z_{jt}}$,译文按原书。)

查看学习笔记对精确 M 步与近似初始坐标的分层推导
查看学习笔记的固定点定位与源文证明审计

最后,$\alpha$ 由分布 $\widehat{\tau}_{\mathrm{marg}}$ 在所有数据点上的经验均值得到,即

$$ \forall k \in [K] : \alpha_k = \frac{1}{nT} \sum_{t=1}^{T} \sum_{i=1}^{n} \widehat{\tau}_{\mathrm{marg}}(t,i,k). $$

6.3.2 时空图上的信念传播(Belief Propagation Using the Space-time Graph)

最大似然估计量通过求解 $\arg\max_{Z} \mathbb{P}\left(A^{1:T} \mid Z\right)$ 找到使似然最大的成员结构,而这里我们转而尝试为每个节点找到使其边际似然最大的块指派。更准确地说,边际似然 $\psi_{k}^{i}(t)$ 是按后验分布节点 $i$ 在时刻 $t$ 属于块 $k$ 的概率,由

$$ \psi_{k}^{i}(t) = \mathbb{P}\left(z_{it} = k \mid A^{1:T}\right). $$

给出。

然后,对每个时刻 $t$,我们把节点 $i$ 指派到块 $\hat{z}_{it}$,使得

$$ \hat{z}_{it} = \operatorname*{arg\,max}_{k \in [K]} \psi_{k}^{i}(t). $$

为计算边际,我们这样建模:给定节点 $i$ 在时刻 $t$ 的每个邻居 $j$ 都发送一个消息 $\psi_{k}^{i \to j}(t)$,它是“若节点 $j$ 不存在,$i$ 属于社区 $k$ 的概率”的一个估计。

由于图是时序的,我们还必须考虑时间演化。更准确地说,在时刻 $t$,每个节点 $i$ 都收到来自其过去与未来副本的消息,分别记为 $\psi^{i(t-1) \to i(t)}$ 与 $\psi^{i(t) \to i(t+1)}$。

空间消息的更新方程为

$$ \begin{array}{l} \displaystyle \psi_{k}^{i \to j}(t) \propto \left( \sum_{\ell} \pi_{k\ell}\,\psi_{\ell}^{i(t-1) \to i(t)} \right) \left( \sum_{\ell} \pi_{\ell k}\,\psi_{\ell}^{i(t+1) \to i(t)} \right) \\ \displaystyle \qquad \times \prod_{\substack{j \,:\, A_{ij}^{t} = 1 \\ j \neq i}} \sum_{\ell} p_{k\ell}\,\psi_{\ell}^{j \to i}(t) \end{array} $$

其中比例关系隐藏了一个施加归一化条件 $\sum_k \psi_{k}^{i \to j} = 1$ 的因子。更新方程忽略了非边,因为在稀疏网络中,非边可以近似为一个全局相互作用。此外,$\psi_{k}^{i \to j}(t)$ 的更新方程不涉及 $j$ 发送给 $i$ 的消息,以避免任何“回声室”效应——否则信息会在 $i$ 与 $j$ 之间以带噪声的方式被放大(更多细节见 Moore, 2017)。以类似方式,时间消息的更新方程由

$$ \psi_{k}^{i(t) \to i(t+1)} \propto \left( \sum_{\ell} \pi_{k\ell}\,\psi_{\ell}^{i(t-1) \to i(t)} \right) \prod_{j \,:\, A_{ij}^{t} = 1} \sum_{\ell} p_{k\ell}\,\psi_{\ell}^{j \to i}(t) $$

给出,且 $\psi_{k}^{i(t-1) \to i(t)}$ 也有类似的表达式。

信念传播(belief propagation)就是随机初始化各消息,然后用更新方程反复更新它们。这通常以异步方式进行:先均匀随机选取一个节点 $i$ 与一个时刻 $t$,对所有 $j$ 与 $k$ 更新 $\psi_{k}^{i \to j}(t)$,以及 $\psi_{k}^{i(t) \to i(t+1)}$ 与 $\psi_{k}^{i(t-1) \to i(t)}$。收敛后,我们用

$$ \begin{array}{rl} & \displaystyle \psi_{k}^{i}(t) \propto \left( \sum_{\ell} \pi_{k\ell}\,\psi_{\ell}^{i(t-1) \to i(t)} \right) \left( \sum_{\ell} \pi_{k\ell}\,\psi_{\ell}^{i(t+1) \to i(t)} \right) \\ & \displaystyle \qquad\quad \times \prod_{j \,:\, A_{ij}^{t} = 1} \sum_{\ell} p_{k\ell}\,\psi_{\ell}^{j \to i}(t). \end{array} $$

计算每个顶点的边际。

本节的最后,我们注意到当 $\pi = rI_K + \frac{1-r}{K}\mathbb{1}_K \mathbb{1}_K^T$ 时,有

$$ \sum_{\ell} \pi_{k\ell}\,\psi_{\ell}^{i(t-1) \to i(t)} = r\,\psi_{k}^{i(t-1) \to i(t)} + \frac{1-r}{K}, $$

$$ \sum_{\ell} \pi_{k\ell}\,\psi_{\ell}^{i(t+1) \to i(t)} = r\,\psi_{k}^{i(t+1) \to i(t)} + \frac{1-r}{K}, $$

这进一步简化了更新方程。此外,在同质模型中(当 $k = \ell$ 时 $p_{k\ell} = p_{\mathrm{in}}$,否则为 $p_{\mathrm{out}}$),有

$$ \sum_{\ell} p_{k\ell}\,\psi_{\ell}^{j \to i}(t) = \lambda\,\psi_{k}^{j \to i}(t) + \frac{1-\lambda}{K}. $$

6.3.3 作为半监督问题的在线推断(Online Inference as a Semi-supervised Problem)

滞后问题(The lagging problem)

在社区成员随时间变化的框架中,直接应用 6.2.3 节推导的时间聚合谱方法做聚类会失败。事实上,随时间变化的社区成员会导致过去相互作用所提供信息的污染。例如,若节点 $i$ 在时刻 $t_1$ 改变其社区指派,那么在寻找它在时刻 $t > t_1$ 的社区成员时,就不应使用节点 $i$ 在前 $t_1$ 个快照中的相互作用。当各层在时间上相关时,这一滞后问题(lagging problem)使情况更加复杂。为避免这一问题,我们提出对节点标签的在线恢复。更具体地:

  • 在时刻 $t = 1$,我们用一个静态社区检测算法输出 $\hat{z}_{\cdot 1} = \left(\hat{z}_{11}, \cdots, \hat{z}_{n1}\right)$,即由第一个快照 $A^1$ 的观测对初始节点标签 $z_{\cdot 1} = (z_{11}, \cdots, z_{n1})$ 的预测;
  • 在时刻 $t > 1$,我们将使用前 $t$ 个快照 $A^1, \ldots, A^t$ 的观测以及此前的预测 $\hat{z}_{\cdot 1}, \cdots, \hat{z}_{\cdot t-1}$。这将被处理为一个半监督学习问题:前一时刻所做的预测 $\hat{z}_{\cdot t-1}$ 被视作时刻 $t$ 真实节点标记 $z_{\cdot t}$ 的带噪声预言机(noisy oracle)。

由马尔可夫结构,时刻 $t > 1$ 的预测归结为仅用时刻 $t-1$ 与 $t$ 的网络以及此前的预测 $\hat{z}_{\cdot t-1}$ 来预测 $z_{\cdot t}$。这可以解释为一个带预言机的噪声半监督问题(见 5.4 节),其中此前的预测 $\hat{z}_{\cdot t-1}$ 扮演了时刻 $t$ 节点标签的预言机信息的角色。这个预言机是带噪声的,因为它带有两类潜在错误。第一,$\hat{z}_{\cdot t-1}$ 不一定恰好等于完美的社区标记 $z_{\cdot t-1}$。第二,由于节点标签随时间变化,$z_{\cdot t-1}$ 并不精确对应于 $z_{\cdot t}$。

Tips:这是本章与第 5 章的回环点:把“上一步的社区预测”当作噪声预言机后,时序在线推断逐时刻化为第 5 章 §5.4 的带噪声半监督聚类问题,下面的命题 6.3 与式 (6.21)–(6.23) 正是这一对应的落实。

6.3.4 具有马尔可夫社区成员的度校正时序 SBM(Degree-corrected Temporal SBM with Markov Community Memberships)

除了 (6.20) 描述的马尔可夫社区结构之外,为简单起见,我们还假设初始标签与转移是均匀的,即

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

换句话说,一个节点以概率 $\eta \in [0,1]$ 保持其标签,并以概率 $1-\eta$ 均匀随机地选择一个标签。

然后,我们假设两个节点 $i$ 与 $j$ 之间的对相互作用是一个马尔可夫过程,只依赖于社区标记与某些度校正参数 $\theta = (\theta_1, \cdots, \theta_N)$。特别地,

$$ \begin{array}{rl} & \displaystyle \mathbb{P}(A \mid z, \theta) = \prod_{1 \leq i < j \leq N} \mathbb{P}\left(A_{ij}^{1} \mid z_{i1}, z_{j1}, \theta_i, \theta_j\right) \\ & \displaystyle \qquad \prod_{t=2}^{T} \mathbb{P}\left(A_{ij}^{t} \mid A_{ij}^{t-1}, z_{it}, z_{jt}, \theta_i, \theta_j\right). \end{array} $$

我们进一步考虑一个同质模型,其中初始分布由

$$ \mathbb{P}\left(A_{ij}^{1} \mid z_{i1}, z_{j1}, \theta_i, \theta_j\right) = \left\{ \begin{array}{ll} \mu^{\theta_i \theta_j}\left(A_{ij}^{1}\right), & \text{若 } z_{i1} = z_{j1}, \\ \nu^{\theta_i \theta_j}\left(A_{ij}^{1}\right), & \text{其他}, \end{array} \right. $$

给出,转移概率由

$$ \mathbb{P}\left(A_{ij}^{t} = b \mid A_{ij}^{t-1} = a, z_{it}, z_{jt}, \theta_i, \theta_j\right) = \left\{ \begin{array}{ll} P_{ab}^{\theta_i \theta_j}, & \text{若 } z_{it} = z_{jt}, \\ Q_{ab}^{\theta_i \theta_j}, & \text{其他}, \end{array} \right. $$

给出。

与 6.2.3 节类似,度校正的初始分布定义为

$$ \mu^{\theta_i \theta_j} = \binom{1 - \theta_i \theta_j \mu_1}{\theta_i \theta_j \mu_1}, \quad \nu^{\theta_i \theta_j} = \binom{1 - \theta_i \theta_j \nu_1}{\theta_i \theta_j \nu_1}, $$

转移概率矩阵由

$$ 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}, \quad 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}, $$

给出,并有假设 $\min_{i,j}\{\theta_i \theta_j \delta\} \le 1$,其中 $\delta = \max\{\mu_1, \nu_1, P_{01}, Q_{01}\}$。我们对度校正参数做归一化,使得对所有 $k$ 成立 $\sum_i \mathbb{1}(z_{i1} = k)\theta_i = \sum_i \mathbb{1}(z_{i1} = k)$。最后,我们假设转移概率与度校正参数不随时间变化,以避免任何参数可辨识性问题(Matias and Miele, 2017)。

在线最大后验估计量(Online Maximum A Posteriori estimator)

下述命题给出上述在线学习问题的 MAP 估计量的表达式。

命题 6.3 在线学习问题的 MAP 估计量

设 $s \in [K]^n$ 是时刻 $t$ 节点标签的一个带噪声预言机,假设它与观测相互作用 $A$ 独立。定义 $s$ 的错误率为 $\rho = \mathbb{P}\left(s_i \neq \hat{z}_{it}\right)$,并假设该错误率对所有节点相同。上述在线学习问题的最大后验(Maximum A Posteriori)估计量定义为

$$ \hat{z}_{\cdot t} = \operatorname*{arg\,max}_{z \in [K]^{n}} \mathbb{P}\left(z \mid A^{t}, A^{t-1}, s\right) $$

且它是任意使下式最大化的标记 $z \in [K]^n$:

$$ \begin{array}{l} \displaystyle \sum_{\substack{i,j \\ z_i = z_j}} \left\{ \ell_{01}^{\theta_i \theta_j}\left(A_{ij}^{t} - A_{ij}^{t-1}A_{ij}^{t}\right) + \ell_{10}^{\theta_i \theta_j}\left(A_{ij}^{t-1} - A_{ij}^{t-1}A_{ij}^{t}\right) + \ell_{11}^{\theta_i \theta_j}A_{ij}^{t-1}A_{ij}^{t} \right. \\ \displaystyle \qquad\quad \left. - \log \frac{Q_{00}^{\theta_i \theta_j}}{P_{00}^{\theta_i \theta_j}} \right\} + 2\lambda \sum_{i=1}^{n} \mathbb{1}\left(z_i = s_i\right), \end{array} $$

其中 $\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}}$,且 $\lambda = \log \frac{1-\rho}{\rho}$。

证明 命题 6.3

由 Bayes 公式,

$$ \mathbb{P}\left(z \mid A^{t}, A^{t-1}, s, \theta\right) \propto \mathbb{P}\left(A^{t} \mid A^{t-1}, z, s, \theta\right)\mathbb{P}\left(z \mid A^{t-1}, s, \theta\right), $$

其中比例符号隐藏了一个不依赖于 $z$ 的项 $\mathbb{P}\left(A^{t} \mid A^{t-1}, s, \theta\right)$。由于 $\mathbb{P}\left(A^{t} \mid A^{t-1}, z, s, \theta\right) = \mathbb{P}\left(A^{t} \mid A^{t-1}, z, \theta\right)$,按照与命题 6.2 的证明相同的做法,对数似然项 $\log \mathbb{P}\left(A^{t} \mid A^{t-1}, z, \theta\right)$ 可以重写为

$$ \begin{array}{l} \displaystyle \frac{1}{2}\sum_{\substack{i,j \\ z_i = z_j}} \left\{ \ell_{01}^{\theta_i \theta_j}\left(A_{ij}^{t} - A_{ij}^{t-1}A_{ij}^{t}\right) + \ell_{10}^{\theta_i \theta_j}\left(A_{ij}^{t-1} - A_{ij}^{t-1}A_{ij}^{t}\right) \right. \\ \displaystyle \qquad\quad \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\}. \end{array} $$

预言机信息等于

$$ \begin{array}{l} \displaystyle \mathbb{P}(z \mid s) = \prod_{i=1}^{n} \frac{\mathbb{P}\left(s_i \mid z_i\right)}{\mathbb{P}\left(s_i\right)}\mathbb{P}\left(z_i\right) \\ \displaystyle \quad = (1-\rho)^{\left| \left\{ i \in [n] \,:\, z_i = s_i \right\} \right|} \rho^{\left| \left\{ i \in [n] \,:\, z_i \neq s_i \right\} \right|} \left(\frac{1}{K}\right)^{n} \\ \displaystyle \quad = \left(\frac{\rho}{1-\rho}\right)^{\left| \left\{ i \in [n] \,:\, z_i \neq s_i \right\} \right|} (1-\rho)^{n}\left(\frac{1}{K}\right)^{n} \end{array} $$

其中用到了节点标签的均匀性。

查看学习笔记的条件性 MAP 推导

MAP 的连续松弛(Continuous relaxation of the MAP)

为简化接下来的推导,本节限于 $K = 2$ 的情形。

记 $A_{\mathrm{pers}}^{t} = A^{t-1} \odot A^{t}$ 为持久边对应的邻接矩阵,$A_{\mathrm{new}} = A^{t} - A_{\mathrm{pers}}^{t}$ 为新形成边对应的邻接矩阵,$A_{\mathrm{old}} = A^{t-1} - A_{\mathrm{pers}}$ 为时刻 $t-1$ 到 $t$ 之间消失边对应的邻接矩阵。那么,利用与 6.2.3 节相同的 Taylor 展开,MAP 估计量可以近似为

$$ \operatorname*{arg\,min}_{z \in \{-1,1\}^{n}} -z^{T}\left( W - \tau \frac{dd^{T}}{2m} \right)z + \lambda(s - z)^{T}(s - z)\tag{6.21} $$

其中 $W^{t} = \alpha_{01}A_{\mathrm{new}}^{t} + \alpha_{10}A_{\mathrm{old}}^{t} + \alpha_{11}A_{\mathrm{pers}}^{t}$,$\alpha_{ab} = \log \frac{P_{ab}}{Q_{ab}}$,$\tau$ 是分辨率参数,$d_i = \sum_{j=1}^{n} W_{ij}^{t}$,且 $m = \frac{1}{2}\sum_{i=1}^{n} d_i$。

这个最小化问题类似于 5.4 节研究的 DC-SBM 中带噪声半监督聚类的问题。我们也可以提出如下连续松弛:

$$ \hat{x} = \operatorname*{arg\,min}_{\substack{x \in \mathbb{R}^{n} \\ x^{T}Dx = 2m}} -x^{T}Mx + \lambda(s - x)^{T}(s - x), $$

其中 $D = \mathrm{diag}(d_1, \cdots, d_n)$,$M = W - \tau\frac{dd^{T}}{2m}$。该松弛的解由模仿 5.4.2 节的推理确定。特别地,记 $D^{-1/2}\left(-M + \lambda I_n\right)D^{-1/2}$ 的特征分解为

$$ D^{-1/2}\left(-M + \lambda I_n\right)D^{-1/2} = Q\Delta Q^{T} $$

其中 $\Delta = \mathrm{diag}(\delta_1, \ldots, \delta_n)$、$QQ^{T} = I_n$,并令 $b = \lambda Q^{T}s$,则 $\hat{x}$ 满足

$$ \left( -M + \lambda I_n - \gamma_{*}D \right)\hat{x} = \lambda s,\tag{6.22} $$

其中 $\gamma_{*}$ 是显式久期方程(secular equation,Gander et al., 1989)

$$ \sum_{i=1}^{n} \left( \frac{b_i}{\delta_i - \gamma} \right)^{2} - 2m = 0\tag{6.23} $$

的最小解。

这导出算法 19。

算法 19 时变社区的在线聚类(online-ssl)

输入:观测图序列 $A^{1:T} = \left(A^1, \ldots, A^T\right)$;社区个数 $K$;静态图聚类算法(记为 $\mathsf{algo}$);参数 $\alpha_{01}, \alpha_{10}, \alpha_{11}$ 与 $\lambda_1, \ldots, \lambda_T$。

输出:节点标记 $Z = \left(z_{it}\right)$。

初始化:计算 $\hat{z}_{\cdot 1} \gets \mathsf{algo}\left(A^1\right)$。

  1. for $t = 2, \ldots, T$ do
  2. 计算 $W = \alpha_{01}A_{\mathrm{new}}^{t} + \alpha_{10}A_{\mathrm{old}}^{t} + \alpha_{11}A_{\mathrm{pers}}^{t}$;
  3. 计算 $M = W - \frac{dd^{T}}{2m}$,其中 $d_i = \sum_{j=1}^{n} W_{ij}$、$m = \frac{1}{2}\sum_{i=1}^{n} d_i$;
  4. 令 $\gamma^{*}$ 为式 (6.23) 的最小解;
  5. 计算 $\hat{x}$ 为式 (6.22) 的解;
  6. 令 $\hat{z}_{\cdot t} = \mathrm{sign}(\hat{x})$。

数值实验(Numerical experiments)

我们在图 6.6 中比较算法 19 与算法 17(带持久边的谱聚类)以及一个对每个快照单独做谱聚类的算法所得的平均精度。特别地,我们观察到,当 $\eta = 1$(即静态社区结构)时,算法 17 如预期一样极为高效:由于它考虑了所有先前的快照,它在这一情形下优于算法 19。相反,当 $\eta \neq 1$ 时,滞后问题出现,算法 17 在若干快照之后最终精度很差。相反,算法 19 在所有快照上都保持很高的精度。

算法 19(online-ssl)与算法 17、逐快照谱聚类的精度对比(a:η = 1,静态社区,算法 17 最优)
(a) $\eta = 1$。
算法 19(online-ssl)与算法 17、逐快照谱聚类的精度对比(b:η = 0.85,出现滞后问题,算法 17 精度下降而算法 19 保持高精度)
(b) $\eta = 0.85$。
图 6.6 算法 19(online-ssl)在 $\alpha_{01} = 1$、$\alpha_{10} = 0$、$\alpha_{11} = 2$ 下的精度,实验在具有 300 个节点、$K = 2$ 个块(均匀先验)且平稳马尔可夫边演化为 $\mu_1 = 0.05$、$\nu_1 = 0.02$、$P_{11} = 0.7$、$Q_{11} = 0.3$ 的时变马尔可夫分块模型上进行。结果对 25 个合成图取平均,误差棒表示标准误。我们与算法 17(加权 SC,$\alpha = 1$、$\beta = 2$)以及一个对每个快照单独做谱聚类的算法进行比较。
不同常数 λ(0.1、0.5、1、1.5、2、5)下算法 19 的精度:λ 在 [0.1, 1] 内性能相近,λ 过大时精度明显下降
图 6.7 算法 19 在 $\alpha_{01} = 1$、$\alpha_{10} = 0$、$\alpha_{22} = 2$ 与不同 $\lambda$ 取值下的精度。模拟在 $n = 300$、$K = 2$、$\mu_1 = 0.05$、$\mu_2 = 0.02$、$P_{11} = 0.7$、$Q_{11} = 0.3$、$\eta = 0.9$ 的时变马尔可夫分块模型上进行。结果对 25 个合成图取平均,误差棒表示标准误。

在图 6.6 中,我们取 $\lambda_t$ 为常数 0.5,而图 6.7 探索了其他可能的取值。我们观察到,当 $\lambda_t$ 取区间 $[0.1, 1]$ 中的常数时,算法 19 输出的性能相近。另一方面,当 $\lambda$ 过大时,算法 19 给予预言机过多的权重,精度变差。实践中,参数 $\lambda_t$ 的选择可以基于数据来优化,例如基于 $\eta$ 或转移矩阵 $P$ 与 $Q$。此外,直观上应让 $\lambda_t$ 随 $t$ 增大,因为可用的时序数据越多,对预言机的置信度越高。我们将此留作未来工作的课题。

进一步阅读(Further Notes)

关于信念传播技术的更全面描述,我们参考 Decelle et al., 2011; Moore, 2017。信念传播由 Ghasemian et al., 2016 引入动态网络,包含边持久性(link persistence)的模型的推广见 Ghasemian, 2019。类似地,Barucca et al., 2018 研究了一个社区成员马尔可夫演化且含边持久性的模型。尽管他们的相互作用设定受限,他们证明了边持久性会增加社区恢复的难度。

最后,一些模型还允许相互作用参数随时间演化(Xu and Hero, 2014; Bhattacharyya and Chatterjee, 2020)。不过需要注意的是,当成员与相互作用核同时随时间变化时,常会出现可辨识性问题(Matias and Miele, 2017)。