“马尔科夫链的粗粒化”的版本间的差异

来自集智百科 - 复杂系统|人工智能|复杂科学|复杂网络|自组织
跳到导航 跳到搜索
第934行: 第934行:
 
===结合复杂网络配分函数与图机器学习的网络流方法===
 
===结合复杂网络配分函数与图机器学习的网络流方法===
  
−
本节介绍的方法主要参考文章<ref name='coarse-graining'>Zhang, Zhang, et al. "Coarse-graining network flow through statistical physics and machine learning." Nature Communications 16.1 (2025): 1605.</ref>,该方法结合了[[复杂网络上的统计物理]]方法与图神经网络方法,从而对[[复杂网络]]进行粗粒化。其核心思想是,粗粒化前后两个网络的[[配分函数]]应该保持不变,将二者的误差作为损失函数,训练一个[[图神经网络]],从而提炼出如何将大网络转化为粗粒化的小网络的映射方法。之所以将复杂网络的配分函数作为一个重要的评价指标,是因为该函数不仅可以体现网络的拓扑结构特征,也能够反映网络的动力学特征,因而反映整体网络功能。实验证明,这套粗粒化网络的方法可以应用到大量的真实网络中,并得到了良好的效果<ref name='coarse-graining'/>。
+
本节介绍的方法主要参考文章<ref name='coarse-graining'>Zhang, Zhang, et al. "Coarse-graining network flow through statistical physics and machine learning." Nature Communications 16.1 (2025): 1605.</ref>,该方法结合了[[复杂网络上的统计物理]]方法与图神经网络方法,从而对[[复杂网络]]进行粗粒化。其核心思想是,粗粒化前后两个网络的[[配分函数]]应该保持不变,将二者的误差作为损失函数,训练一个[[图神经网络]],从而提炼出如何将大网络转化为粗粒化的小网络的映射方法。之所以将[[复杂网络的配分函数]]作为一个重要的评价指标,是因为该函数不仅可以体现网络的拓扑结构特征,也能够反映网络在多个时间尺度上的动力学特征,因而反映整体网络功能。实验证明,这套粗粒化网络的方法可以应用到大量的真实网络中,并得到了良好的效果<ref name='coarse-graining'/>。
  
 
下面介绍具体的方法。首先,我们用[[扩散方程]] <math>\frac{d\psi}{d\tau}=-L\psi</math> 来描述网络的扩散过程,其中[[拉普拉斯矩阵]] <math>L=D-A</math> 描述网络的拓扑结构,<math>A</math> 是邻接矩阵,<math>D</math> 是[[度矩阵]]。扩散方程的解为 <math>\psi_{\tau}=e^{-\tau L}\psi_0</math>,体现了信息随时间的演化,其中 <math>e^{-\tau L}</math> 为时间演化算子,<math>[e^{-\tau L}]_{ij}</math> 表示经过时间 <math>\tau</math> 后从节点 <math>j</math> 到节点 <math>i</math> 的流量。
 
下面介绍具体的方法。首先,我们用[[扩散方程]] <math>\frac{d\psi}{d\tau}=-L\psi</math> 来描述网络的扩散过程,其中[[拉普拉斯矩阵]] <math>L=D-A</math> 描述网络的拓扑结构,<math>A</math> 是邻接矩阵,<math>D</math> 是[[度矩阵]]。扩散方程的解为 <math>\psi_{\tau}=e^{-\tau L}\psi_0</math>,体现了信息随时间的演化,其中 <math>e^{-\tau L}</math> 为时间演化算子,<math>[e^{-\tau L}]_{ij}</math> 表示经过时间 <math>\tau</math> 后从节点 <math>j</math> 到节点 <math>i</math> 的流量。

2026年7月14日 (二) 09:06的版本

马尔科夫链是一类特殊的随机过程:在给定当前状态的信息后,对未来状态的预测不依赖更早的过去状态。对于一个包括 [math]\displaystyle{ N }[/math] 个状态的 马尔科夫链,所有状态之间的概率转移构成一个转移概率矩阵 [math]\displaystyle{ P }[/math]。该矩阵有 [math]\displaystyle{ N }[/math] 行 [math]\displaystyle{ N }[/math] 列,且满足概率归一化条件(所有元素位于 0 到 1 之间,所有行的和为 1)。本词条限定在 平稳马尔科夫链 上讨论,也就是 [math]P[/math] 不随时间和状态而变。

所谓马尔科夫链的粗粒化,是指找到一个新的马尔科夫链,它具有 [math]\displaystyle{ 1\lt M\lt N }[/math] 个状态,其对应一个 [math]\displaystyle{ M }[/math] 行 [math]\displaystyle{ M }[/math] 列概率转移矩阵 [math]\displaystyle{ P' }[/math]。它同样满足概率归一化条件,可以看作对 [math]\displaystyle{ P }[/math] 所描述动力学过程的一种近似的抽象、聚合和约简。

文献中着重强调对马尔科夫链做粗粒化的意义包括:

  1. 当我们研究一个状态极多的复杂系统时,往往不能关注每一个微观状态。在粗粒化中,我们希望过滤掉微观尺度中不感兴趣的噪声和异质性,从中总结出一些中尺度或宏观尺度的规律;
  2. 有些状态的转移概率非常相似,所以可以被看成同一类状态。对这些状态进行划分和归并可以减少马尔科夫链的冗余性,得到一个更简洁有效的表示;
  3. 在用到马尔科夫决策过程的强化学习算法中,对马尔科夫链做粗粒化可以减少状态空间的大小,提高训练过程的样本效率和任务表现[1]。

对任意马尔科夫链的转移概率矩阵做约简,存在无穷多种可能的做法。然而在实际粗粒化过程中,我们必须从这些方法中做筛选,使约简后的矩阵 [math]P'[/math] 尽可能保留原矩阵的主要特征。为刻画这些特征是否被保留,文献中提出了不同的标准。总体原则包括:

  1. 首先,我们想尽量降低动力学在压缩后的维度;
  2. 其次,我们想在压缩过程中尽量保留动力学中的信息;
  3. 最后,我们希望压缩后的过程和原本的过程尽量保持一致。

在这个词条中,我们会围绕着这三个原则展开讨论,并介绍一系列针对马尔科夫链的粗粒化方法。


粗粒化

马尔科夫链不只是一个马尔科夫矩阵。它还涉及不同的状态,以及状态之间的转移概率,所以我们可以从不同角度来看待粗粒化:

  1. 对马尔科夫矩阵做粗粒化:在看待一个方阵时,我们可以把它想象成一个网络的邻接矩阵。节点为状态,连边为状态转移概率,状态转移的过程为一个粒子在节点间跳转,而粗粒化则是对这个网络做降维。
  2. 对状态空间做粗粒化:我们这里讨论的是离散状态的马尔科夫链,状态空间的大小等于马尔科夫链的大小。粗粒化可以看作对状态空间做分组并投影到新的状态空间。
  3. 对概率空间做粗粒化:从概率论出发,把状态的发生看作事件,通过计算事件发生的频率来构建概率空间。在降维表达事件的同时也构建了新的概率空间。

这三种出发点彼此相关,但具体关注的对象不同:第一种关注状态之间的连边拓扑结构,第二种关注不同状态的历史分布和相似性,第三种关注如何降维表达整个系统中的事件概率分布。


对马尔科夫链做粗粒化并没有统一定义,目前该主题也没有完整的文献综述。在许多文献中,粗粒化(coarse-graining)、聚类/聚合(partitioning/lumping/clustering/aggregation)和降维(dimension reduction)往往是重叠乃至近似等价的。

本词条只讨论离散时间马尔科夫链,且动力学 [math]\displaystyle{ P }[/math] 是时齐(不随时间变化)的,不讨论连续马尔科夫过程。我们定义微观状态空间为[math]\displaystyle{ S=\{s_1, s_2, ... ,s_n\} }[/math],微观转移函数[math]\displaystyle{ P \in \mathcal{R}_{n \times n} }[/math]是一个大小为 n 的马尔科夫矩阵。同样,粗粒化后的宏观状态空间为[math]\displaystyle{ S'=\{s'_1, s'_2, ... ,s'_r\} }[/math],宏观转移函数[math]\displaystyle{ P' \in \mathcal{R}_{r \times r} }[/math]。多数文献中会出现如图1的定义[2]:

图1:马尔科夫链一步状态转移过程与粗粒化过程示意图

对状态空间 [math]\displaystyle{ S }[/math] 和转移函数 [math]\displaystyle{ P(t):S \rightarrow S }[/math],我们需要一个状态分组的投影函数 [math]\displaystyle{ \Phi:S \rightarrow S' }[/math],和一个在新状态空间上的转移函数 [math]\displaystyle{ P'(t):S' \rightarrow S' }[/math],而且 [math]\displaystyle{ S' }[/math] 和[math]\displaystyle{ P' }[/math] 要符合马尔科夫链定义。

如果把这个降维/粗粒化过程拆开看,它可以分为三步。

  1. 划分(对状态空间进行粗粒化):对状态进行分组,得出[math]\displaystyle{ \Phi }[/math],决定 [math]\displaystyle{ S }[/math] 空间中的状态要如何被投影到 [math]\displaystyle{ S' }[/math] 空间上。
  2. 聚合(对动力学/转移概率进行粗粒化):基于分组函数 [math]\displaystyle{ \Phi }[/math],我们要决定如何把原来的状态转移矩阵 [math]\displaystyle{ P }[/math] 降维成新的 [math]\displaystyle{ P' }[/math]。
  3. 判断粗粒化方案 [math]\displaystyle{ \Phi }[/math] 和 [math]\displaystyle{ P' }[/math] 符合什么样的性质,从而决定是否接受当前的划分和聚合方案。


划分:对状态空间进行粗粒化

我们假设投影函数 [math]\displaystyle{ \Phi }[/math] 是一个 [math]\displaystyle{ n \times r }[/math] 的线性投影矩阵,且符合下列两个条件:(1)所有元素非负;(2)每行加和等于1。

马尔科夫链的粗粒化可以分成硬划分(hard partitioning)和软划分(soft partitioning)。硬划分是指把若干个微观状态分成一个宏观状态类,且一个微观状态不能同时属于多个宏观状态类;软划分则可能出现一个微观状态同时属于多个宏观类的情况。

在硬划分的情况下,[math]\displaystyle{ \Phi_{i, j} = \begin{cases} 1\text{, if}\ s_i \in s_j', \\ 0\text{, otherwise} \end{cases} }[/math],即一个表示 [math]\displaystyle{ i }[/math] 属于 [math]\displaystyle{ j }[/math] 类的零一分组矩阵。如果是软划分,[math]\displaystyle{ \Phi }[/math] 则不一定是零一矩阵[3]。

在这个词条里,我们主要讨论硬划分[4]。对软划分感兴趣的朋友可以自行参阅这两篇论文[5] [6]。


聚合:对转移概率进行粗粒化

得出 [math]\displaystyle{ \Phi }[/math] 后,对微观马尔科夫矩阵 [math]\displaystyle{ P }[/math] 的降维可以表示为

[math]\displaystyle{ \begin{aligned} P' = W P \Phi, \end{aligned} }[/math]

 

 

 

 

(1)

其中 [math]\displaystyle{ W \in \mathbb{R}^{r \times n} }[/math],其中 [math]\displaystyle{ W }[/math] 的任意元素满足:[math]\displaystyle{ W_{ij}\geq 0, \forall i,j\in\{1,2,\cdots,n\} }[/math],并要求转换后的 [math]\displaystyle{ P' }[/math] 满足行归一化条件。Buchholz[3]的文章中把 [math]\displaystyle{ W }[/math] 称作分配器矩阵(distributor matrix),而 [math]\displaystyle{ \Phi \in \mathbb{R}^{n \times r} }[/math] 称作收集器矩阵(collector matrix)。

我们可以这样理解这个公式:宏观状态从 [math]\displaystyle{ t }[/math] 到 [math]\displaystyle{ t+1 }[/math] 时刻的过程 [math]\displaystyle{ P' }[/math],可以表达为从图1左上角的 [math]\displaystyle{ t }[/math] 时刻宏观状态出发,先经过一个“逆向粗粒化”操作 [math]\displaystyle{ W }[/math] 得到 [math]\displaystyle{ t }[/math] 时刻微观状态,再经微观动力学 [math]\displaystyle{ P }[/math] 得到 [math]\displaystyle{ t+1 }[/math] 时刻的微观状态,最后经粗粒化 [math]\displaystyle{ \Phi }[/math] 得到 [math]\displaystyle{ t+1 }[/math] 时刻的宏观状态。

从式子1可以看到,[math]\displaystyle{ \Phi }[/math] 和 [math]\displaystyle{ P }[/math] 都是给定的,所以这一步要做的是决定 [math]\displaystyle{ W }[/math]。[math]\displaystyle{ W }[/math] 的基本限制是使 [math]\displaystyle{ P' }[/math] 符合马尔科夫矩阵定义;在此基础上,[math]\displaystyle{ W }[/math] 仍有无限多种选项,我们可以从中挑出更符合目标特征的方案。


粗粒化方案的要求和优化目标

目前文献中尚没有一个严格、完备的马尔科夫链粗粒化定义。在此,本词条将简要梳理对马尔科夫链粗粒化操作提出的一些要求,这些基本要求可分为约束条件和优化指标两个方面。

问题拆分

粗粒化并不只是“把若干状态聚成几类”的聚类问题。我们需要同时考虑多个层面:保留什么尺度上的信息、如何保证宏观动力学仍然自洽、以及如何在压缩和保真之间做权衡。我们需要区分:哪些微观差异只是会被快速混合掉的局部噪声、快变量或冗余通路,哪些差异则决定了系统未来的主要演化方向、稳态概率流以及跨尺度的因果结构。

更直观地说,一个好的宏观态应当把“对外行为几乎一样”的微观状态并在一起,把“长期命运明显不同”的状态分开。后文提到的成块性、因果态、尺度闭合、有效信息和动力学可逆性,分别对应这个问题的不同侧面:有的在回答“哪些状态适合合并”,有的在回答“合并后是否仍然自洽”,有的则在回答“哪个尺度最有信息价值”。

下图展示了我们对粗粒化问题的拆分:

图2:对粗粒化问题的拆分


粗粒化与时间演化算子的可交换性

首先,一个重要的粗粒化标准是:粗粒化操作和时间演化操作应当是可交换(commutative)的两个算符。以图1所示,我们从某个微观状态 [math]\displaystyle{ s^{(t)} }[/math] 出发,到达下一刻宏观状态 [math]\displaystyle{ A^{(t+1)} }[/math] 有两条不同路径[2]:

  1. 先经过微观动力学的时间演化达到 [math]\displaystyle{ s^{(t+1)} }[/math],再进行粗粒化得到 [math]\displaystyle{ A^{(t+1)} }[/math];
  2. 先通过状态的粗粒化得到 [math]\displaystyle{ A^{(t)} }[/math],再经过宏观动力学的演化到达 [math]\displaystyle{ A^{(t+1)} }[/math]。

如果我们把粗粒化和(宏观/微观)动力学看作是两个算子,那这两个算子的可交换性指的是:从任意一个微观状态分布 [math]\displaystyle{ s^{(t)} }[/math] 出发,走这两条路径得到的 [math]\displaystyle{ A^{(t+1)} }[/math] 的分布是一样的,即

[math]\displaystyle{ \begin{aligned} s^{(t)} \Phi P' = s^{(t)} P \Phi \end{aligned} }[/math]

 

 

 

 

(2)

这两条路径的一致性在粗粒化领域中非常重要,也在成块性和因果抽象等工作中被反复提及 [7] [2] [3] [8] [9] [10] [11] 。满足一致性的粗粒化方案能够充分保留动力学中的有效信息,并丢弃无用信息。具体来说,宏观状态已经捕捉了与未来演化相关的所有“慢变量”或“相关信息”,而被粗粒化掉的那些微观细节,只影响瞬时噪声或快速弛豫过程,不会长期、系统性地改变宏观演化。

对于马尔科夫链,粗粒化时我们通常希望得到一个新的马尔科夫过程:宏观态的转移只依赖当前宏观态,而不依赖更久远的历史。式子(2)说明,如果粗粒化满足可交换性,那么宏观过程可以写成一个真正意义上的马尔科夫链:拥有一个明确的宏观转移矩阵 [math]\displaystyle{ P' }[/math],而且宏观态的演化不需要记忆“刚才具体是哪个微观态导致了这个宏观态”。如果不满足,那宏观态的演化一般会带有记忆效应,未来的转移概率还会依赖这个宏观态是如何形成的,也就是依赖过去的轨迹。


成块性

上一节说的可交换性给出了粗粒化最理想的要求:无论先在微观演化再粗粒化,还是先粗粒化再在宏观演化,最后得到的宏观分布都应当一致。但这个要求本身还不是一个方便检查的划分条件。对于马尔科夫链来说,成块性(Lumpability)正是可交换性的一个具体判据。我们在马尔科夫链的成块性词条中已证明,成块性和“对任意初始分布成立的精确可交换粗粒化”是等价的。

马尔科夫链的成块性最早由 Kemeny 和 Snell 在其 1969 年著作 Finite Markov Chains[7] 中阐述。Lump 一词有“块”的意思。判断一条马尔科夫链是否可成块,就是判断是否存在一种状态划分,能够把原始链压缩成一条保留块间转移机制的简化链。换句话说,成块性回答了,哪些微观状态可以被合并,而不会让宏观过程产生额外的记忆效应。这种划分需确保在抹除聚类内部状态的微观差异后,仍能精确保持块间转移机制。因此,成块性为状态空间的合理划分提供了严格的数学判据,也成为马尔科夫链维度约简中最常用的标准之一。

成块性的原始定义比较复杂,详见马尔科夫链的成块性。这里我们介绍其充分必要条件,作为一个更直观的解析。在一些文献中,作者会直接使用这一通俗版解释,而且它与原定义之间存在充分必要条件关系[7]。

成块性的充分必要条件

设 [math]\displaystyle{ p_{s_k \rightarrow s_m} = p(s^{(t)} = s_m | s^{(t-1)} = s_k) }[/math] 为微观状态从 [math]t-1[/math] 时刻 [math]\displaystyle{ s_k }[/math] 到 [math]t[/math] 时刻 [math]\displaystyle{ s_m }[/math] 的转移(条件)概率;

设 [math]\displaystyle{ p_{s_k \rightarrow A_i} = p(s^{(t)} \in A_i | s^{(t-1)} = s_k) = \sum_{s_m \in A_i} p_{s_k \rightarrow s_m} }[/math] 为 [math]t-1[/math] 时刻的微观状态 [math]\displaystyle{ s_k }[/math] 到 [math]t[/math] 时刻的宏观状态 [math]\displaystyle{ A_i }[/math] 的转移(条件)概率(定义为从[math]\displaystyle{ s_k }[/math]到 [math]\displaystyle{ A_i }[/math] 内所有微观状态的条件概率之和);

设 [math]\displaystyle{ p_{A_i \rightarrow A_j} = p(A^{(t)} = A_j | A^{(t-1)} = A_i) }[/math] 为宏观状态从 [math]t-1[/math] 时刻 [math]\displaystyle{ A_i }[/math] 到 [math]t[/math] 时刻 [math]\displaystyle{ A_j }[/math] 的转移(条件)概率。

作者[7]提出了判断一个马尔科夫链对给定划分 [math]\displaystyle{ A }[/math] 是否成块的充分必要条件为:

对于任意一对 [math]\displaystyle{ A_i, A_j }[/math],每一个属于 [math]\displaystyle{ A_i }[/math] 的状态 [math]\displaystyle{ s_k }[/math] 的 [math]\displaystyle{ p_{s_k \rightarrow A_j} }[/math] 都相同。也就是说,[math]\displaystyle{ A_i }[/math] 里的每一个状态 [math]\displaystyle{ s_k }[/math] 转移到 [math]\displaystyle{ A_j }[/math] 中所有状态的转移概率之和是相同的。

[math]\displaystyle{ \forall i, j, \forall s_k, s_m \in A_i, p_{s_k \rightarrow A_j} = p_{s_m \rightarrow A_j} }[/math]

 

 

 

 

(3)

图3:成块性充分必要条件实例

满足这个条件后,同一宏观组内的微观状态对所有宏观组的总转移概率都相同,因此可以把这些微观状态合并为一个宏观状态。合并后,宏观状态之间的转移概率可以由组内任意一个微观状态到目标宏观组的总转移概率给出。

例如,在上图中,[math]\displaystyle{ s_1 }[/math]和[math]\displaystyle{ s_2 }[/math]都属于[math]\displaystyle{ A_1 }[/math],它们到[math]\displaystyle{ A_1 }[/math]和[math]\displaystyle{ A_2 }[/math]的转移概率相同;同样,[math]\displaystyle{ s_3 }[/math]和[math]\displaystyle{ s_4 }[/math]都属于[math]\displaystyle{ A_2 }[/math],它们到各宏观状态的转移概率也相同。因此,[math]\displaystyle{ [[1,2], [3, 4]] }[/math]是[math]\displaystyle{ P }[/math]的一个成块划分。

基于这个条件,我们可以定义成块划分的宏观动力学[math]\displaystyle{ A^{(t)} = A^{(t-1)} P' }[/math]:

对于给定的微观状态 [math]\displaystyle{ S }[/math],微观动力学 [math]\displaystyle{ p_{s_k \rightarrow s_m} }[/math] 和成块的划分 [math]\displaystyle{ A }[/math],宏观动力学 [math]\displaystyle{ p_{A_i \rightarrow A_j} }[/math] 等于 [math]\displaystyle{ A_i }[/math] 中任意状态 [math]\displaystyle{ s_k }[/math] 到[math]\displaystyle{ A_j }[/math] 的概率:

[math]\displaystyle{ \begin{aligned} p_{A_i \rightarrow A_j} = p_{s_k \rightarrow A_j} = \sum_{s_m \in A_j} p_{s_k \rightarrow s_m}, k \in A_i \end{aligned} }[/math]

 

 

 

 

(4)

这个公式表明,群组 [math]\displaystyle{ A_i }[/math] 到群组 [math]\displaystyle{ A_j }[/math] 的转移概率 = 群组 [math]\displaystyle{ A_i }[/math] 中任意状态 [math]\displaystyle{ s_k }[/math] 到群组 [math]\displaystyle{ A_j }[/math] 的转移概率 = 群组 [math]\displaystyle{ A_i }[/math] 中任意状态 [math]\displaystyle{ s_k }[/math] 到群组 [math]\displaystyle{ A_j }[/math] 中的所有状态的转移概率之和。

我们也可以直接定义 [math]\displaystyle{ P }[/math] 和 [math]\displaystyle{ P' }[/math] 之间的关系:

[math]\displaystyle{ P' = \Phi^T P (\Phi^T)^{-1} = \Phi^T P \Phi (\Phi^T \Phi)^{-1} }[/math],

其中,

[math]\displaystyle{ \Phi^T \Phi=\textrm{diag}(\{\sum_{k=1}^n\mathbb{1}_{s_k\in A_i}\}_{1\leq i\leq n}) }[/math]

是一个对角矩阵,对角线上的元素为每个 [math]\displaystyle{ A_i }[/math] 包括了多少个微观状态的计数,其逆矩阵就是把对角线上的元素全部取倒数。


在讨论成块性时,必须明确状态划分的方式。不同的划分方式可能导致不同的结论。一个马尔科夫链可能在某种划分下是成块的,而在另一种划分下则不成块。比如上图例子中,分组 [math]\displaystyle{ \{\{s_1,s_2\}, \{s_3,s_4\}\} }[/math] 是成块的,但是当我们考虑分组 [math]\displaystyle{ \{\{s_1,s_2,s_3\}, \{s_4\}\} }[/math] 时,组内状态的 [math]\displaystyle{ p_{s_k \rightarrow A_j} }[/math] 并不相同,因此不满足成块性的条件定义。

成块性的主要性质

成块性作为一种较为严格的划分方式约束,虽然在压缩维度上可能不如其他方法显著,但它能够保留在宏观层面上所有有用的信息,并且确保粗粒化前后过程的一致性。

宏观最大预测性

宏观最大预测性是指,在将微观状态粗粒化为宏观状态后,我们只用宏观状态就能预测后续的宏观状态。

具体来说,只需要知道宏观状态 [math]\displaystyle{ A^{(t)} }[/math],就可以最大化 [math]\displaystyle{ A^{(t+1)} }[/math] 的预测效果,而微观状态 [math]\displaystyle{ s^{(t+1)} }[/math] 虽然能提供一些额外信息,但是对于预测 [math]\displaystyle{ A^{(t+1)} }[/math] 并无帮助。

这意味着条件互信息 [math]\displaystyle{ I(A^{(t+1)}; s^{(t+1)} | A^{(t)}) = 0 }[/math]。宏观动力学不需要任何来自微观状态的修正,就能独立且正确地反映系统在宏观层面上的演化。

尽管成块性划分能够准确预测宏观动力学,但它无法确保像计算力学中的因果态那样的微观最大预测性。因为在聚类过程中会损失微观状态的信息,使得我们无法辨别某个[math]\displaystyle{ A^{(t)} }[/math]具体是哪个微观态[math]\displaystyle{ s^{(t)} }[/math]。这个信息的损失是不可逆的,所以从微观粗粒化到宏观之后,就无法再还原到微观了。

成块划分和粗粒化的可交换性

如本节开头所述,在强成块的意义下,成块性与对任意初始分布成立的精确可交换粗粒化等价。因此,成块划分不仅给出一种状态聚合方式,也保证聚合后的宏观动力学可以独立闭合。详细推导可见马尔科夫链的成块性词条。

一致性

我们在马尔科夫链的成块性中推导了结论:基于成块性划分的粗粒化方法构建的宏观网络与微观网络是一致的。

所谓的一致性,请参考:动力学的一致性检验词条


因果态

上一节的成块性关心的是:在给定状态划分后,宏观过程是否仍然是一条封闭的马尔科夫链。它要求同一宏观块内的微观状态,对各个宏观块具有相同的总转移概率,因此保证的是宏观层面的可预测性:知道当前宏观态,就足以预测下一步宏观态。

因果态则是一个更强的预测性概念。它并不只要求同组状态对各个宏观块的转移概率相同,而是要求这些状态对整个未来(微观)状态序列具有相同的预测分布。也就是说,因果态把所有“对未来预测等价”的历史归为一类。对于一阶马尔科夫链,如果两个微观状态属于同一个因果态,那么它们的未来分布相同,因而自然也会满足成块性所要求的宏观转移一致;但反过来,一个成块划分未必是因果态划分,因为它可能只保留了宏观转移概率,而丢掉了预测微观未来所需的信息。所以,因果态是较成块性更严格的概念。

此部分为对因果态的简略介绍,因果态相关详情请参照计算力学词条。

定义

在计算力学中,状态被理解为智能体接收到的全部历史信息,因此系统从过去到未来的变化被描述为一个离散的平稳随机过程。在一般情形下,过去的历史都可能影响对未来的预测。因此,为保证通用性,因果态原来的定义适用于无限长的历史序列上,而不预设系统具有马尔科夫性。

本文讨论的是马尔科夫链的状态粗粒化。若系统是一个 [math]\displaystyle{ K }[/math] 阶马尔科夫链([math]\displaystyle{ t }[/math] 时刻之后的状态分布由前 [math]\displaystyle{ K }[/math] 个时刻的状态决定),则预测未来只需保留最近 [math]\displaystyle{ K }[/math] 个状态,可以把这段历史片段看作新的状态。对于本文的一阶马尔科夫链,这段历史片段就是当前微观状态 [math]\displaystyle{ s_t }[/math]。记所有可能的微观状态为 [math]\displaystyle{ S }[/math],下一时刻状态为 [math]\displaystyle{ s_{t+1}\in S }[/math],因果态就可以直接定义在 [math]\displaystyle{ S }[/math] 上。

因果态这种宏观态划分方法,具有最强的预测能力。形式化定义如下:

对于任意微观状态 [math]\displaystyle{ s_i }[/math] 和 [math]\displaystyle{ s_j }[/math],若给定当前状态 [math]\displaystyle{ s_i }[/math] 和 [math]\displaystyle{ s_j }[/math] 时,下一时刻状态的分布完全相同,则称 [math]\displaystyle{ s_i }[/math] 与 [math]\displaystyle{ s_j }[/math] 等价,也叫预测等价性(predictive equivalence),记为:

[math]\displaystyle{ s_i \sim s_j \triangleq P(S_{t+1}=s_k|S_t=s_i ) = P(S_{t+1}=s_k | S_t=s_j ),\quad \mathrm{for~all~}s_k\in S }[/math]

 

 

 

 

(4)


计算力学认为,包括人类在内的智能体,每时每刻都在对周围环境观察,得到大量的观测数据,并在自己的内部(大脑中)构建简洁的机制来预测观测数据的变化(一种隐马尔科夫模型)。所以,引入了 ϵ-machine(ϵ-机器)的概念来描述这种机制。ϵ-机器根据当前时间步的状态,来预测下一时刻状态。其工作机制是:根据 [math]\displaystyle{ s_t }[/math] 而判断 [math]\displaystyle{ t }[/math] 时刻系统处于哪个因果态 [math]\displaystyle{ E_t }[/math] ,并根据该因果态预测下一时刻状态 [math]\displaystyle{ s_{t+1} }[/math]。

因果态可定义为ϵ-机器的映射 [math]\displaystyle{ \epsilon{:}S\mapsto2^S }[/math],[math]\displaystyle{ E = \{E_t\}_{t \in \mathbb{N} } }[/math],[math]\displaystyle{ E_t = \epsilon( s_t ) }[/math],满足:[math]\displaystyle{ \epsilon(s_i)\equiv\{s_j|\mathrm{P}(S_{t+1}=s_k\mid S_t=s_i)=\mathrm{P}(S_{t+1}=s_k\mid S_t=s_j),\mathrm{for~all~}s_k\in S,s_j\in S\} }[/math]

所有的因果态的集合由预测等价关系"[math]\sim[/math]"诱导的商集定义:[math]\mathcal{S}=S/\sim[/math]


因果态的主要性质

为什么我们要关注因果态这样一种特定的宏观态?下面是因果态的两个主要性质,它们体现了因果态为何是一种理想的宏观态:

微观最大预测性

微观最大预测性是指,在将微观状态粗粒化为宏观状态后,我们只用宏观状态就能预测后续的微观状态。

具体来说,只需要知道宏观状态 [math]\displaystyle{ A^{(t)} }[/math],就可以最大化 [math]\displaystyle{ s^{(t+1)} }[/math] 的预测效果。

对于宏观态 [math]\displaystyle{ \mathcal{\forall R} }[/math],都有 [math]\displaystyle{ H[S_{t+1}|\mathcal{R}]\geq H[S_{t+1}|\mathcal{S}] }[/math],[math]\displaystyle{ S_{t+1} }[/math] 为下一时刻的微观状态变量,[math]\displaystyle{ \mathcal{S} }[/math] 为因果态划分。[math]\displaystyle{ H[S_{t+1}|\mathcal{R}] }[/math] 和 [math]\displaystyle{ H[S_{t+1}|\mathcal{S}] }[/math] 是 [math]\displaystyle{ S_{t+1} }[/math] 的条件熵。可以理解为因果态集合 [math]\displaystyle{ \mathcal{S} }[/math] 是在所有的划分宏观态集合 [math]\displaystyle{ \mathcal{R} }[/math]中,预测能力最强的一种宏观态,证明过程请参照计算力学词条。

最小随机性

对于所有的宏观态 [math]\displaystyle{ \hat{\mathcal{R}} }[/math] 和 [math]\displaystyle{ \hat{\mathcal{R}}^{\prime} }[/math],都有 [math]\displaystyle{ H[\hat{\mathcal{R}}^{\prime}|\hat{\mathcal{R}}]\geq H[\mathcal{S}^{\prime}|\mathcal{S}] }[/math],其中 [math]\displaystyle{ \hat{\mathcal{R}}^{\prime} }[/math] 和 [math]\displaystyle{ \mathcal{S}^{\prime} }[/math] 分别是该过程的下一时刻状态和下一时刻因果态。

该性质可以被理解为在相同预测能力的前提下,因果态集合 [math]\displaystyle{ \mathcal{S} }[/math] 在所有的宏观态 [math]\displaystyle{ \mathcal{R} }[/math] 的集合中,随机性最小。

用互信息的角度去理解的话,上式等价于 [math]\displaystyle{ I(\mathcal{S}^{\prime};\mathcal{S})\geq I(\hat{\mathcal{R}}^{\prime};\hat{\mathcal{R}}) }[/math],可以理解为因果态为所有的与它自己的下一时刻的互信息中的最大的一个宏观态。


尺度闭合性

尺度闭合性讨论的是:粗粒化后的宏观动力学,能否在不反复回到底层微观细节的情况下,自洽地演化、预测和解释。换句话说,我们关心的是:在给定粗粒化 [math]\displaystyle{ Z_t=f(X_t) }[/math] 后,宏观尺度是否已经“足够独立”,可以单独承担预测、干预和机制描述的任务[12]。关于尺度闭合性的详细讨论请见不同尺度之间的闭合性。

跨尺度因果态和[math]\displaystyle{ \nu }[/math]-机器

前文已经介绍过因果态和 ϵ-机器。在引入闭合性概念之前,再介绍一种机器,[math]\nu[/math]-机器,它专门处理“微观历史预测宏观未来”。

设微观过程为 [math]\displaystyle{ X_t }[/math],粗粒化后的宏观过程为 [math]\displaystyle{ Z_t=f(X_t) }[/math]。宏观的 [math]\epsilon[/math]-机器是用宏观历史 [math]\displaystyle{ \overleftarrow{Z_t} }[/math] 去预测宏观未来 [math]\displaystyle{ \overrightarrow{Z_t} }[/math];而 Rosas 等人提出的[math]\nu[/math]-机器,则是用微观历史 [math]\displaystyle{ \overleftarrow{X_t} }[/math] 去预测宏观未来 [math]\displaystyle{ \overrightarrow{Z_t} }[/math]。它在微观历史上定义的等价关系为:

[math]\displaystyle{ \overleftarrow{x_t} \equiv \overleftarrow{x_{t'}} \ \ \text{iff}\ \ p(\overrightarrow{z_t}|\overleftarrow{x_t}) = p(\overrightarrow{z_{t'}}|\overleftarrow{x_{t'}}) }[/math]

 

 

 

 

(5)

也就是说,若两段微观历史对宏观未来给出的预测分布完全相同,它们就属于同一个跨尺度因果态。把这些等价类记作 [math]\displaystyle{ U_t=υ(\overleftarrow{X_t}) }[/math]。直观地看,[math]\nu[/math]-机器保留的是“对宏观未来真正有用的那部分微观信息”,因此可以把它理解为一种面向跨尺度预测的最优压缩。

在给定粗粒化 [math]\displaystyle{ Z_t=f(X_t) }[/math] 后,可以同时出现三种相关状态表征:

  1. 微观 [math]\epsilon[/math]-机器的因果态 [math]\displaystyle{ E_t }[/math],保留了对微观未来最完整的预测信息;
  2. 宏观 [math]\epsilon[/math]-机器的因果态 [math]\displaystyle{ E'_t }[/math],只保留宏观历史中的信息;
  3. 跨尺度 [math]\nu[/math]-机器的因果态 [math]\displaystyle{ U_t }[/math],保留了微观历史中对宏观未来仍然有用的信息。

三种闭合

因果闭合

若宏观过程 [math]\displaystyle{ Z_t }[/math] 的 [math]\epsilon[/math]-机器因果态 [math]\displaystyle{ E'_t }[/math] 与跨尺度 [math]\nu[/math]-机器的因果态 [math]\displaystyle{ U_t }[/math] 之间存在一一对应关系,那么称该粗粒化满足因果闭合(Causal Closure)。

它的物理意义是:微观层面上那些会影响宏观未来的因果差异,已经被宏观历史完整表达出来了。因此,在宏观层面实施干预,与在微观层面做出对应干预后再投影到宏观层,二者在宏观结果上是等价的。若因果闭合成立,我们就不必为了分析宏观干预效果而不断回到底层微观细节。

信息闭合

信息闭合(Information Closure)强调的不是“干预”,而是“预测”。它要求宏观历史与微观历史对宏观未来提供的信息量相同,即

[math]\displaystyle{ I(\overleftarrow{Z_t};\overrightarrow{Z_t}) = I(\overleftarrow{X_t};\overrightarrow{Z_t}) }[/math]

 

 

 

 

(6)

这意味着:若我们已经知道宏观历史 [math]\displaystyle{ \overleftarrow{Z_t} }[/math],再额外补充微观历史 [math]\displaystyle{ \overleftarrow{X_t} }[/math],也不会提升对宏观未来的预测能力。它的物理意义是,宏观尺度在信息上已经“自给自足”,微观细节对宏观预测来说是冗余的。

Anil Seth等人在文章[13]中提出了“动力学解耦”(Dynamical independence)的概念,并认为符合信息闭合的系统中的宏观和微观动力学已经实现了“解耦”。一个人在跑步的时候,其身体内的细胞会进行各种各样的运动和反应。然而,我们不需要知晓所有这些细胞的状态,仅凭他的肢体动作就能判断他下一秒是否会继续跑步。也就是说,细胞的信息对我们来说是冗余的,而细胞这个微观层面和肢体这个宏观层面的动力学则实现了“解耦”。

对给定的粗粒化函数 [math]\displaystyle{ f }[/math] 而言,因果闭合和信息闭合是等价的[12]。前者从干预与因果的角度表达同一件事,后者从预测与信息的角度表达同一件事。

计算闭合

计算闭合更关注“机制能否一起被粗粒化”。若存在一个对应于机制层的粗粒化函数 [math]\displaystyle{ f^* }[/math],使得微观 ϵ-机器的因果态 [math]\displaystyle{ E_t }[/math] 可以被粗粒化成宏观 ϵ-机器的因果态 [math]\displaystyle{ E'_t }[/math],即 [math]\displaystyle{ E'_t=f^*(E_t) }[/math],那么称该粗粒化满足计算闭合(Computational Closure)。

它的物理意义是:宏观动力学不只是状态层面上的简化,也可以看作微观计算机制的一种简化版。也就是说,先找微观机制再粗粒化,与先粗粒化再找宏观机制,在原则上是兼容的。


最后,三种闭合性彼此之间的关系可以概括为:对同一个粗粒化 [math]\displaystyle{ Z_t=f(X_t) }[/math],因果闭合 = 信息闭合,而它们又都会推出计算闭合;反过来则未必成立。因此,计算闭合是三者中相对更宽松的条件。

尺度闭合 与 成块性

图4:强成块性 < 信息闭合/因果闭合 < 计算闭合;强成块性 < 弱成块性

成块性和尺度闭合都关心“粗粒化后是否仍然自洽”,但它们关注的角度不同。成块性是马尔科夫链的内部结构的条件,强调宏观过程是否仍然是一个真正的马尔科夫链;尺度闭合则更强调跨尺度的预测、干预和机制是否能够独立成立。

需要先区分两种成块性。本文前面默认讨论的是强成块性:给定划分后,无论微观初始分布如何,只要两个微观状态属于同一个宏观块,它们到任意宏观块的总转移概率都相同。因此,强成块性保证宏观过程对任意初始分布都是封闭的一阶马尔科夫链。弱成块性则更宽松:它不要求所有初始分布都成立,只要求在某些特定初始分布(常见的是稳态分布或与划分兼容的一类分布)下,粗粒化后的宏观过程仍然表现为马尔科夫性。因此,弱成块性可以包含一些不满足强成块条件、但在特定观测分布下仍能形成宏观马尔科夫动力学的情形。

Geiger 和 Temmel 的工作[14],以及 Rosas 等人的工作[12],给出了成块性与尺度闭合性之间的联系。粗略地说,强成块性是最干净、最严格的情形:它不仅保证宏观马尔科夫性,也通常落在信息闭合/因果闭合之内。信息闭合和因果闭合关心的是“微观历史是否还能为宏观未来提供额外信息”,因此比强成块性覆盖更广;计算闭合又更宽松,只要求机制层面的粗粒化能够对应起来。弱成块性与尺度闭合不是简单的包含关系。一方面,弱成块性比强成块性宽松,可能与信息闭合/因果闭合有重叠;另一方面,它依赖特定初始分布或特定观测条件,因此不能推出所有尺度闭合性质。

可以把二者的关系理解为:成块性更像“单尺度的马尔科夫动力学闭合”,尺度闭合性则更像“跨尺度的预测、干预和机制自洽”。强成块性要求最强,结论最干净;弱成块性和尺度闭合则覆盖更多中间情形,更适合讨论实际系统中的近似涌现、跨尺度预测和机制简化。关于该概念更详细丰富的讨论请见不同尺度之间的闭合性。


因果涌现:最大化有效信息

因果涌现(causal emergence)是指动力系统中的一类特殊涌现现象,即系统在宏观尺度展现出更强的因果特性。特别地,对于一类马尔科夫动力系统来说,在对其状态空间进行适当粗粒化以后,形成的宏观动力学会展现出比微观更强的因果特性,也就是更大的有效信息,那么我们称该系统发生了因果涌现[15][16]。

定义

在马尔科夫链中,任意时刻的状态变量[math]X_t[/math]都可以看作是原因,而下一时刻的状态变量[math]X_{t+1}[/math]就可以看作是结果,这样马尔科夫链的状态转移矩阵就是它的因果机制。

有效信息(Effective Information,简称EI)是因果涌现 Causal Emergence理论中的一个核心概念,它可以用来衡量一个马尔科夫动力学的因果效应的强度。这里,因果效应是指把动力学看做一个黑箱,那么不同的输入分布就会导致不同的输出分布,二者之间联系的紧密程度就是因果效应:

[math]\displaystyle{ EI = I(X_t,X_{t+1}|do(X_t\sim U(\mathcal{X})))=\frac{1}{N}\sum_{i=1}^N D_{KL}(P_i||\bar{P}) }[/math]

 

 

 

 

(EI)

其中 [math]\displaystyle{ \bar{P} }[/math] 是转移矩阵各行的平均转移向量。这个表达式说明,EI 本质上度量的是:不同初始状态所对应的转移分布,与“平均转移分布”之间到底有多可区分。若各行彼此差异很大,动力学的因果效应就更强;若各行非常相似,EI 就会较小。

下表展示的是三个不同的转移概率矩阵,以及它们的EI数值:

马尔科夫链示例
[math]\displaystyle{ P_1=\begin{pmatrix} &0 &0 &1 &0& \\ &1 &0 &0 &0& \\ &0 &0 &0 &1& \\ &0 &1 &0 &0& \\ \end{pmatrix} }[/math],

[math]\displaystyle{ P_2=\begin{pmatrix} &1/3 &1/3 &1/3 &0& \\ &1/3 &1/3 &1/3 &0& \\ &0 &0 &0 &1& \\ &0 &0 &0 &1& \\ \end{pmatrix} }[/math],

[math]\displaystyle{ P_3=\begin{pmatrix} &1/4 &1/4 &1/4 &1/4& \\ &1/4 &1/4 &1/4 &1/4& \\ &1/4 &1/4 &1/4 &1/4& \\ &1/4 &1/4 &1/4 &1/4& \\ \end{pmatrix} }[/math].

[math]EI(P_1)=2[/math] bits [math]EI(P_2)=1[/math] bits [math]EI(P_3)=0[/math] bits

上面所列 [math]\displaystyle{ P_1, P_2, P_3 }[/math] 的EI各为:2比特、1比特和0比特。由此可见,如果 [math]\displaystyle{ P_1, P_2, P_3 }[/math] 中出现更多的独热(one-hot)向量,则EI值就会更大。也就是说,如果在状态转移的过程中,从某一时刻到下一时刻的跳转越确定,则EI值就会倾向于越高。



最大化有效信息的划分

前面提到,对马尔科夫链做粗粒化的一个动机,就是想抹除一些微观动力学上的噪音,从而得到在宏观上更具规律性和因果性的动力学。所以,我们可以利用 EI 这一马尔科夫链的因果效应度量,来衡量聚类后宏观动力学的因果强度。

但是,EI 的大小和状态空间大小有关,这一性质在比较不同尺度的马尔科夫链时很不方便。我们需要一个尽可能不受尺度效应影响的因果效应度量,因此需要对有效信息 EI 做归一化处理,得到一个和系统尺寸无关的量化指标。根据Erik Hoel和Tononi等人的工作,应使用均匀分布即最大熵分布下的熵值 [math]\displaystyle{ \log N }[/math] 作为分母对 EI 进行归一化,这里的[math]N[/math]为状态空间[math]\mathcal{X}[/math]中的状态数[15]。那么归一化后的 EI 为:

[math]\displaystyle{ Eff=\frac{EI}{\log N} }[/math]

进一步定义归一化指标也称为有效性(effectiveness)。


有了有效信息这一度量指标后,我们便可以讨论马尔科夫链上的因果涌现。对于一个马尔科夫链,观察者可以建立多尺度视角来区分微观和宏观。首先,原始的马尔科夫概率转移矩阵 [math]\displaystyle{ P_N }[/math] 定义了微观动力学;其次,经过对微观态的粗粒化映射(coarse-graining)后,观察者可以得到对应的宏观状态(即分组),这些宏观状态彼此之间的概率转移可以用宏观的马尔科夫概率转移矩阵 [math]\displaystyle{ P_M }[/math] 来刻画。对两个动力学分别计算 EI 后,如果宏观的 EI 大于微观 EI,我们就说该系统发生了因果涌现。

图5:因果涌现示意图

可以定义一个新的指标直接度量因果涌现的程度:

[math]\displaystyle{ ce = Eff(P_M) - Eff(P_N) = \frac{EI(P_M)}{\log M} - \frac{EI(P_N)}{\log N} }[/math]

这里[math]P[/math]为微观状态的马尔科夫概率转移矩阵,维度为:[math]N\times N[/math],这里N为微观的状态数;而[math]P'[/math]为对[math]P[/math]做粗粒化操作之后得到的宏观态的马尔科夫概率转移矩阵,维度为[math]M\times M[/math],其中[math]M<N[/math]为宏观状态数。如果计算得出的[math]\displaystyle{ ce \gt 0 }[/math],则称该系统发生了因果涌现,否则没有发生。

下面,我们展示一个具体的因果涌现的例子:

马尔科夫链示例

[math]\displaystyle{ P_m=\begin{pmatrix} &1/3 &1/3 &1/3 &0& \\ &1/3 &1/3 &1/3 &0& \\ &1/3 &1/3 &1/3 &0& \\ &0 &0 &0 &1& \\ \end{pmatrix} }[/math],

[math]\displaystyle{ P_M=\begin{pmatrix} &1 &0 & \\ &0 &1 & \\ \end{pmatrix} }[/math].

[math]\begin{aligned}&Det(P_m)=0.81\ bits,\\&Deg(P_m)=0\ bits,\\&EI(P_m)=0.81\ bits\end{aligned}[/math] [math]\begin{aligned}&Det(P_M)=1\ bits,\\&Deg(P_M)=0\ bits,\\&EI(P_M)=1\ bits\end{aligned}[/math]

在这个例子中,微观态的转移矩阵是一个 4×4 的矩阵,其中前三个状态彼此以 1/3 的概率相互转移。这导致该转移矩阵的确定性较小,因此 EI 也不高,仅为 0.81。对该矩阵进行粗粒化时,我们把前三个状态合并为一个状态 a,把最后一个状态视为另一个宏观态 b。这样,原本三个微观态之间的转移就变成了宏观态 a 的内部转移,因此转移概率矩阵变成了[math]P_M[/math],其 EI 为 1。在这个例子中,可以计算其因果涌现度量为:

[math]\displaystyle{ ce=Eff(P_M)-Eff(P_m)=1-0.405=0.595 }[/math]

当[math]\displaystyle{ ce }[/math]越大,说明宏观动力学的因果效应更强,也说明宏观动力学能够更好地表示系统的因果涌现。

动力学一致性

Erik Hoel 和 Klein 在其工作[17]中提出了动力学一致性,衡量了粒子分别在原始马尔科夫链和粗粒化后链中演化得到的状态轨迹的一致程度。动力学一致性在某种程度上反映了动力学与粗粒化算子的可交换性,但比可交换性更灵活、更宽松,可以计算“不可交换但近似交换”的程度,因此可以作为粗粒化方案合理性的一种判断依据。

以下关于动力学的一致性检验的描述节选自复杂网络中的因果涌现词条:

定义微观网络 [math]\displaystyle{ G }[/math] 与宏观网络 [math]\displaystyle{ G_M }[/math],分别在两个网络上进行随机游走。在某个时间 [math]\displaystyle{ t }[/math],[math]\displaystyle{ G }[/math] 上的预期分布为 [math]\displaystyle{ P_m(t) }[/math],而 [math]\displaystyle{ G_M }[/math] 上的预期分布为 [math]\displaystyle{ P_M(t) }[/math]。

将 [math]\displaystyle{ P_m(t) }[/math] 分布叠加到宏观网络 [math]\displaystyle{ G_M }[/math] 的相同节点上,得到 [math]\displaystyle{ P_{M|m}(t) }[/math] 分布。通过计算 [math]\displaystyle{ P_M(t) }[/math] 和 [math]\displaystyle{ P_{M|m}(t) }[/math] 之间的 KL 散度来衡量其不一致性(inconsistency)。如果KL散度结果为零,则表示动力学一致。公式如下:

[math]\displaystyle{ inconsistency=\sum_{t=0}^T D_{KL}[P_M(t)||P_{M|m}(t)] }[/math]

也就是说,不一致性表示宏观和微观网络节点在不同时刻概率分布的KL散度之和。

近似动力学可逆性与SVD

所谓的近似动力学可逆性定义为概率转移矩阵(TPM)的奇异值之和([math]\alpha=1[/math]),因此基于近似动力学可逆性的粗粒化方法就是利用TPM的SVD分解,截取前k个大奇异值,并往这k个大奇异值所对应的奇异向量所组成的子空间做投影的方法。因此,这种粗粒化方法的原则就是尽量保证粗粒化后系统的近似动力学可逆性保持不变。

动力学可逆性

所谓的严格意义下的动力学可逆性(Dynamical Reversibility),是指一个马尔可夫链的转移矩阵[math]P[/math]本身是矩阵可逆的,且逆矩阵满足行归一化特性,可以证明,当[math]P[/math]是一个置换矩阵(Permutation Matrix)的时候,它严格地满足动力学可逆性的要求,但是对于一般的马尔可夫转移矩阵则不满足。因此,我们寻求刻画一个TPM接近一个严格可逆置换阵的方法。

近似动力学可逆性

为了刻画一个动力学TPM[math]P[/math]是否与一个置换矩阵或动力学可逆矩阵足够接近,文献[18]提出了近似动力学可逆性(Approximate Dynamical Reversibility,记为:[math]\displaystyle{ \Gamma_{\alpha} }[/math])的概念,用于量化接近可逆动力学的程度。

定义(近似动力学可逆性):对于给定的TPM:[math]P[/math],它的 [math]\displaystyle{ \alpha }[/math] 阶近似动力学可逆性定义为:

[math]\displaystyle{ \begin{aligned} \Gamma_{\alpha}=\sum_{i=1}^{N}\sigma_{i}^{\alpha}\end{aligned} }[/math]

其中 [math]\displaystyle{ \sigma_i }[/math] 是 [math]\displaystyle{ P }[/math] 的奇异值。[math]\displaystyle{ \Gamma_{\alpha} }[/math] 越大,说明动力学越接近可逆、冗余通路越少。当 [math]\displaystyle{ \alpha \to 0 }[/math] 时,[math]\displaystyle{ \Gamma_0 }[/math] 衡量的是 [math]\displaystyle{ P }[/math] 的(非)简并性;而 [math]\displaystyle{ \alpha \to 2 }[/math] 时,[math]\displaystyle{ \Gamma_2 }[/math] 衡量的是 [math]\displaystyle{ P }[/math] 的确定性。通过调整 [math]\displaystyle{ \alpha \in (0, 2) }[/math] 来调整考虑确定性和简并性对最终可逆性度量的影响。更进一步的范数性质和技术证明,请详见基于可逆性的因果涌现理论。

基于SVD的粗粒化方法

如果将TPM视作一个通信信道,则近似动力学可逆性 [math]\Gamma_{\alpha}[/math] 反映了该通讯信道的信息流通的能力,TPM中的每一个行向量对应为一个信息通路(information pathway)。理想化的 [math]P[/math] 的粗粒化方案应该能尽可能地保留这种能力,也就是说即使在减少了该通讯信道的有效通道数的前提下,仍然能够保留信道整体的近似可逆性程度。由于 [math]\Gamma_{\alpha}=\sum_{i=1}^N\sigma_i^{\alpha}[/math],其中 [math]\sigma_i[/math] 为第 [math]i[/math] 大奇异值,因此,该粗粒化方案可以通过保留前 [math]k[/math] 个奇异值的方式来作为核心原则(这里 [math]k[/math] 为宏观TPM中保留下来的状态数)。

文献[18]中给出了一种实施该方案的具体算法,算法的基本流程如下:

  1. 给定含有 [math]\displaystyle{ N }[/math] 个状态的马尔科夫矩阵 [math]\displaystyle{ P }[/math],对其做 SVD 分解 [math]\displaystyle{ P = U \Sigma V^T }[/math];
  2. 根据奇异值谱设定阈值 [math]\displaystyle{ \epsilon }[/math],保留前 [math]\displaystyle{ r_\epsilon }[/math] 个右奇异向量,构成 [math]\displaystyle{ V_{N \times r_\epsilon} }[/math];
  3. 对 [math]\displaystyle{ P }[/math] 的行向量做投影:[math]\displaystyle{ \tilde{P} = P \cdot V_{N \times r_\epsilon} }[/math],得到每个状态在“核心通路”子空间中的表示;
  4. 对 [math]\displaystyle{ \tilde{P} }[/math] 的行向量用 K-Means 聚类为 [math]\displaystyle{ r_\epsilon }[/math] 组,得到划分 [math]\displaystyle{ \Phi }[/math];
  5. 基于划分 [math]\displaystyle{ \Phi }[/math] 和稳态流矩阵 [math]\displaystyle{ F = \mathrm{diag}(\pi) P }[/math],通过 [math]\displaystyle{ F' = \Phi^T F \Phi }[/math] 构建宏观稳态流矩阵,再对每行归一化得到宏观转移矩阵 [math]\displaystyle{ P' }[/math]。

这里的投影 [math]\displaystyle{ \tilde{P} }[/math] 可以理解为:只保留转移矩阵中最主要的几条信息通道,把冗余或近乎无效的方向滤掉,再按照这些主导模式来归组状态。最后用 [math]\displaystyle{ F' = \Phi^T F \Phi }[/math] 而不是直接对 [math]\displaystyle{ P }[/math] 做行列聚合,是为了保证粗粒化与时间演化尽可能可交换,从而使宏观动力学仍然是合法的马尔科夫链。

与因果涌现理论的关系

事实上,基于动力学可逆性和SVD的粗粒化方法最一开始是作为一种替代最大化有效信息的因果涌现量化方法,即基于可逆性的因果涌现理论而提出的,它是一种量化因果涌现的新框架。

相似点

事实上,之所以该理论可以作为对最大化有效信息因果涌现方案的近似替代,是因为对于任意的时齐马尔可夫过程,近似动力学可逆性的对数近似等价于有效信息([math]\log\Gamma_{\alpha}\sim EI[/math]),即指标[math]\Gamma_{\alpha}[/math]与有效信息高度相关,因此它也可以用来刻画动力学的因果效应强度,因此因果涌现其实也可以近似等价于近似动力学可逆性的涌现。

我们可以通过最大化EI来寻找最优的马尔可夫链的粗粒化方案,自然也可以通过最大化近似动力学可逆性[math]\Gamma_{\alpha}[/math]来寻找最优的粗粒化方案,只不过在动力学可逆性理论中,这个最大化求解过程可以直接通过对转移概率矩阵进行奇异值分解而得到。

基于近似动力学可逆性和SVD这套粗粒化方案的合理性,来源于其与“最大化有效信息”的因果涌现定义之间的密切关系。在Erik Hoel 因果涌现理论的框架内,给定马尔科夫动力学系统的因果涌现大小依赖于粗粒化策略的选择,而最优粗粒化策略可以通过最大化宏观动力学的 EI 来寻找。

可以证明,要获得最大化有效信息的粗粒化策略,近似等价于采用向最大奇异值对应的奇异向量方向投影的策略。严格推导可见基于有效信息的因果涌现理论。由于对离散马尔科夫链来说,粗粒化策略往往是由比较严格的 0 或 1 构成的向量,因此不能简单地按奇异向量直接投影。但是,如果粗粒化策略尽可能与最大的几个奇异向量平行,就可以让粗粒化后动力学的 EI 尽可能更大。因此,最大化 EI 的必要条件近似可以理解为:向最大奇异值对应的奇异向量空间做投影。在实际应用中,粗粒化策略 [math]\displaystyle{ \phi }[/math] 的约束(例如宏观动力学 TPM 的分组和归一化)会阻止 [math]\displaystyle{ \phi_j }[/math] 与奇异向量完美对齐。因此,基于 SVD 的粗粒化策略只是最大化 EI 的近似必要条件。

如果把一个马尔科夫动力学看作通信信道,TPM中的每一行看作一个信息通路,那么系统发生因果涌现就等价于其中存在冗余的信息通路,也就是某些状态转移彼此相似或线性相关。消除这些冗余通路的最好办法,就是利用 SVD 找到等于零或接近于零的奇异值及其对应的奇异向量,这实际上正是最大化有效信息的精髓所在。事实上,可以证明,最大化宏观动力学有效信息的粗粒化方案,近似等价于将系统动力学投影到较大奇异值对应的奇异向量方向上的方案[18]。而系统所能达到的最大信息传输效率提升(归一化的近似动力学可逆性),就可以定义为因果涌现的强度。

差异性

基于可逆性的因果涌现理论与其他因果涌现理论最大的区别在于:它不需要指定粗粒化策略,仅从马尔科夫转移矩阵的奇异值谱就能判断因果涌现是否发生,并直接定义清晰因果涌现(clear causal emergence)和模糊因果涌现(vague causal emergence)的概念。这在一定程度上规避了因果涌现结果依赖粗粒化策略、以及最大化有效信息难以优化求解的问题。其次,这套理论指出,一个动力系统的因果涌现特性仅由系统线性化算符的奇异值谱决定,因而是一套更能反映系统动力学特性本身的因果涌现理论。相关详情请参照基于可逆性的因果涌现理论词条。

根据式EI,EI定义为[math]P[/math]中所有行向量彼此之间的相似性程度(以KL散度来衡量)。因而,通过EI最大化得到的粗粒化方案可以通过消除这些彼此相似的行向量,从而实现概率转移矩阵的粗粒化约简。但是仅考虑这种用KL散度衡量的行向量相似性是不够的,下面两个 4×4 马尔科夫链展示了另一种行向量上的冗余性:

[math]\displaystyle{ P_{\mathrm{lin}}=\begin{pmatrix} 0.9 & 0.1 & 0 & 0 \\ 0 & 0 & 0.1 & 0.9 \\ 0.27& 0.03& 0.07& 0.63\\ 0.63& 0.07& 0.03& 0.27 \end{pmatrix} }[/math]

[math]\displaystyle{ P_{\mathrm{cmp}}=\begin{pmatrix} 0.9 & 0.1 & 0 & 0 \\ 0 & 0 & 0.1 & 0.9 \\ 0.63& 0.03& 0.07& 0.27\\ 0.27& 0.07& 0.03& 0.63 \end{pmatrix} }[/math]

在 [math]\displaystyle{ P_{\mathrm{lin}} }[/math] 中,第三行和第四行并不是新的独立通路,而是前两行的线性组合:

[math]\displaystyle{ P_3=0.3P_1+0.7P_2,\qquad P_4=0.7P_1+0.3P_2 }[/math]

其物理含义是,当当量随机游走粒子以30%概率位于第一个状态,以70%概率位于第2个状态的时候,它们的概率转移就与处于第三状态的随机游走粒子没有任何区别了,因此第三状态下的概率转移通路与以30%走第一状态下的概率转移通路和以70%概率走第二状态的概率转移通路没有区别,存在着冗余。同样的道理,对于第四行也有类似的冗余性。

而 [math]\displaystyle{ P_{\mathrm{cmp}} }[/math] 使用了同样的行分布形状和同样的平均输出分布,但第三、四行不再只是前两行的线性混合。换句话说,两个矩阵在 EI 看重的分布差异上几乎被设计成一样,但在行向量彼此之间的线性相关性上则不同。诸如[math]P_{\mathrm{lin}}[/math]这样的线性相关性引起的相似性用EI就无法刻画,但却可以被近似动力学可逆性指标[math]\Gamma_{\alpha}[/math]来刻画。

所以,EI无法刻画这种由线性相关带来的冗余性,但是近似动力学可逆性则可以。

对上面两个矩阵进行计算,可以得到:

EI 与动力学可逆性对比
矩阵 EI 非零奇异值 [math]\displaystyle{ \Gamma_0 }[/math] [math]\displaystyle{ \Gamma_1 }[/math] [math]\displaystyle{ \Gamma_2 }[/math]
[math]\displaystyle{ P_{\mathrm{lin}} }[/math] [math]\displaystyle{ 0.559 }[/math] bits [math]\displaystyle{ 1.281,\ 0.975 }[/math] [math]\displaystyle{ 2 }[/math] [math]\displaystyle{ 2.256 }[/math] [math]\displaystyle{ 2.591 }[/math]
[math]\displaystyle{ P_{\mathrm{cmp}} }[/math] [math]\displaystyle{ 0.559 }[/math] bits [math]\displaystyle{ 1.281,\ 0.972,\ 0.074 }[/math] [math]\displaystyle{ 3 }[/math] [math]\displaystyle{ 2.327 }[/math] [math]\displaystyle{ 2.591 }[/math]

这两个矩阵的 EI 完全相同,因为它们的平均输出分布相同,各行转移分布的熵也相同。也就是说,从 EI 的角度看,它们具有相同的因果效应强度。但从奇异值谱看,[math]\displaystyle{ P_{\mathrm{lin}} }[/math] 只有两条独立的信息通路,而 [math]\displaystyle{ P_{\mathrm{cmp}} }[/math] 有三条独立的信息通路。因此,[math]\displaystyle{ \Gamma_0 }[/math] 和 [math]\displaystyle{ \Gamma_1 }[/math] 能区分这两种动力学结构,而 EI 无法区分。这个例子说明,[math]\displaystyle{ \Gamma_{\alpha} }[/math] 的优势不在于替代 EI,而在于它能刻画 EI 不敏感的线性冗余:有些状态转移在概率分布上并不相同,却可能只是已有信息通道的线性混合。需要注意的是,[math]\displaystyle{ \Gamma_2 }[/math] 在这个例子中相同,因为它更接近于衡量矩阵的确定性或能量大小;真正揭示通路冗余的是更偏向秩结构的 [math]\displaystyle{ \Gamma_0 }[/math],以及介于二者之间的 [math]\displaystyle{ \Gamma_1 }[/math]。

稳态分布下的随机动力学划分

随机动力学为我们提供了另一个视角来讨论马尔科夫链的 状态划分,这一状态划分是以状态分布为基础的。

当我们观察马尔科夫链的一个具体样本时,其在每一步都处于某个确定状态。而所谓状态分布,是指给定初始分布 [math]\displaystyle{ \mu }[/math] 后,系统在 [math]\displaystyle{ t }[/math] 时刻落在各个状态的概率。若存在分布 [math]\displaystyle{ \pi }[/math] 满足 [math]\displaystyle{ \pi P = \pi }[/math],则称[math]\pi[/math]为稳态分布。所以,当系统在到达稳态后,再继续演化,整体分布会保持不变。

我们可以把马尔科夫链放进随机动力系统(Random Dynamical System)的框架中来理解(详细的数学公式和证明可参考马尔科夫链中的随机动力学)。在这个框架里,不变性(Invariance)是理解稳态分布的关键。稳态分布描述的是分布层面的不变,而随机吸引子描述的则是状态集合层面的不变。一个状态集合若要成为随机吸引子,通常需要满足三点:

  1. 紧致性 (Compactness):对有限状态马尔科夫链,这一点天然满足;
  2. 不变性 (Invariance):系统从集合中的状态出发演化后,仍落在下一时刻的吸引子中;
  3. 吸引性 (Pullback Attraction):从足够久远的过去出发,系统最终会被拉回到该集合附近。

因此,随机吸引子刻画的是系统长期反复停留、且不会轻易逃离的结构。正是在这个意义上,我们可以在稳态下讨论可约、周期、滑过态这些概念和分组。

基础划分

从长期行为看,马尔科夫链中的状态可以先分成两类:滑过态(Transient)和常返态(Recurrent)。滑过态只是在过渡阶段被访问,离开后最终可能不再返回;常返态则会被反复访问。常返态进一步组成闭合常返类:同一类中的状态彼此可达,而且一旦进入该类就不会再离开。若一个常返类的返回步数存在公因子 [math]\displaystyle{ d\gt 1 }[/math],它还可以继续分解为 [math]\displaystyle{ d }[/math] 个周期子类,系统按固定顺序在这些子类之间循环。

随机吸引子和对应分组

在稳态或长时间视角下,随机吸引子给出的不是任意聚类,而是一种按照长期命运进行的自然划分:会进入同一个闭合常返类的状态、处在同一个周期相位中的状态,或者只承担过渡作用的滑过态,应该被区别对待。若这些集合被分作不同的宏观态,我们就可以尝试构造宏观转移矩阵。

在稳态视角下,马尔科夫链常见有三种情形:

  • 不可约、非周期:整个状态空间形成唯一的平凡吸引子 [math]\displaystyle{ A(\omega)=X }[/math]。这说明系统的长期活动遍及整个空间,因此随机动力学本身并不给出更细的自然分组。比如

[math]\displaystyle{ P=\begin{pmatrix} 1/2&1/2&0\\ 1/4&1/2&1/4\\ 0&1/2&1/2 \end{pmatrix} }[/math],

这个链不可约且有自环,所以是非周期的。随机吸引子是全集 [math]\displaystyle{ X=\{1,2,3\} }[/math]。若按吸引子粗粒化,只有一个宏观态 [math]\displaystyle{ A=X }[/math],对应的宏观转移矩阵为 [math]\displaystyle{ P'=(1) }[/math]。这是一种平凡粗粒化,说明不可约非周期链本身没有由吸引子诱导出的非平凡分组。

  • 可约:设 [math]\displaystyle{ S=T \cup R_1 \cup \cdots \cup R_m }[/math],其中 [math]\displaystyle{ T }[/math] 是滑过态集,[math]\displaystyle{ R_i }[/math] 是闭合常返类。此时每个 [math]\displaystyle{ R_i }[/math] 都对应一个确定性的吸引子 [math]\displaystyle{ A_i=R_i }[/math]。系统一旦进入某个闭合常返类,就不会再离开。比如

[math]\displaystyle{ P=\begin{pmatrix} 1/3&1/3&1/3&0&0\\ 1/3&1/3&1/3&0&0\\ 1/3&1/3&1/3&0&0\\ 0&0&0&1&0\\ 0.6&0&0&0.4&0 \end{pmatrix} }[/math],

其中 [math]\displaystyle{ R_1=\{1,2,3\} }[/math] 和 [math]\displaystyle{ R_2=\{4\} }[/math] 是两个闭合常返类,[math]\displaystyle{ T=\{5\} }[/math] 是滑过态。若为了描述过渡阶段而保留滑过态,按 [math]\displaystyle{ \{T,R_1,R_2\} }[/math] 的顺序粗粒化,可以得到

[math]\displaystyle{ P'=\begin{pmatrix} 0&0.6&0.4\\ 0&1&0\\ 0&0&1 \end{pmatrix} }[/math]。

这个例子中,宏观态 [math]\displaystyle{ R_1 }[/math] 和 [math]\displaystyle{ R_2 }[/math] 都是吸引子:一旦进入就出不去;滑过态 [math]\displaystyle{ T }[/math] 只描述进入吸引子之前的过渡。若只看稳态分布,[math]\displaystyle{ T }[/math] 的概率质量最终为 0,系统只会停留在某个闭合常返类内;此时稳态流归一化只能在具有正稳态质量的闭合常返类上定义,例如在 [math]\displaystyle{ \{R_1,R_2\} }[/math] 上得到吸收型宏观动力学。

  • 不可约、周期为 [math]\displaystyle{ d\gt 1 }[/math]:状态空间可分解为 [math]\displaystyle{ X=C_0 \cup C_1 \cup \cdots \cup C_{d-1} }[/math],并按 [math]\displaystyle{ C_0 \to C_1 \to \cdots \to C_{d-1} \to C_0 }[/math] 循环。此时吸引子不再是固定不动的闭合集合,而是随时间在周期子类之间切换的非平凡随机吸引子。比如

[math]\displaystyle{ P=\begin{pmatrix} 0&0&1/2&1/2\\ 0&0&1/2&1/2\\ 1/2&1/2&0&0\\ 1/2&1/2&0&0 \end{pmatrix} }[/math] ,

该链不可约,但周期为 2。状态空间可分为 [math]\displaystyle{ C_0=\{1,2\} }[/math] 和 [math]\displaystyle{ C_1=\{3,4\} }[/math],系统每一步都在两个周期类之间来回跳转。按 [math]\displaystyle{ \{C_0,C_1\} }[/math] 粗粒化后,宏观转移矩阵为 [math]\displaystyle{ P'=\begin{pmatrix} 0&1\\ 1&0 \end{pmatrix} }[/math]。

因此,周期情形下的“吸引子分组”不是“进入某个吸引子后永远停在里面”,而是保留系统的周期相位:从一个周期类必然转到下一个周期类。若只看稳态分布,两个周期类的总概率可以保持不变;但从单条轨道看,宏观态仍然会周期性跳转。

所以,常返类、周期类和滑过态之间的边界提供了一种值得考虑的粗粒化基准,按照此边界分组可以把长期命运相同的状态显式区分出来。

基于稳态流的聚合方法

虽然随机吸引子划分可以看作一种粗粒化,但无法保证其划分是成块的。上述的三个例子中的分组均不是成块的。所以,像成块性那样直接平均 [math]\displaystyle{ P_{ij} }[/math] 等于强行给组内(可能不均匀的)节点分配同样的权重,会导致微观和宏观动力学的不一致。

更自然的做法是先使用稳态流矩阵 [math]\displaystyle{ F_{ij}=\pi_iP_{ij} }[/math] 聚合概率流:[math]\displaystyle{ F'_{ab}=\sum_{i\in A_a}\sum_{j\in A_b}\pi_iP_{ij} }[/math]。此时 [math]\displaystyle{ F'_{ab} }[/math] 是从宏观态 [math]\displaystyle{ A_a }[/math] 流向 [math]\displaystyle{ A_b }[/math] 的稳态概率流量,还不是转移概率。若 [math]\displaystyle{ \pi'_a=\sum_{i\in A_a}\pi_i\gt 0 }[/math],则对应的宏观转移概率应为 [math]\displaystyle{ P'_{ab}=\frac{F'_{ab}}{\pi'_a} =\frac{\sum_{i\in A_a}\sum_{j\in A_b}\pi_iP_{ij}}{\sum_{i\in A_a}\pi_i} }[/math]。这样相当于按照状态在稳态时的出现频率安排权重,可以保证该宏观动力学在稳态时和微观动力学保持一致。

滑过态在稳态分布下概率质量为 0,即 [math]\displaystyle{ \pi'_a=0 }[/math],无法通过稳态流归一化定义其宏观出流。因此,滑过态通常更适合作为过渡区单独处理,或者按照其到各吸引子的吸收概率和转移模式进一步细分。详细的数学公式和证明可参考马尔科夫链中的随机动力学。


不同粗粒化方法之间的异同

对上述的概念进行一个总结:

  1. 粗粒化与时间演化算子的可交换性
    • 本质:要求 / 指标(既是理想要求,又可度量偏差)
  2. 成块性
    • 本质:一种结构性性质
    • 作用:给出在“完全可交换、完全宏观马尔科夫”的极限下,粗粒化应满足的严格条件;也可以当作一种“硬约束目标”。
  3. 因果态
    • 本质:结构 + 性质,“单尺度下实现最大预测性的最优表征”的理论工具。
  4. 尺度闭合性(及因果/信息/计算闭合)
    • 本质:要求 / 性质,是从“多尺度因果结构”角度,对一个粗粒化是否“自洽”的刻画。
  5. 因果涌现:最大化有效信息
    • 本质:优化目标 + 评价指标(EI)
  6. 近似动力学可逆性
    • 本质:优化目标:平均每通道近似动力学可逆性([math]\Gamma_{\alpha}/N[/math])+ SVD分解
  7. 稳态分布下的随机动力学划分
    • 本质:一种自然划分方式 + 结构性质,像是从随机动力系统视角导出的一类“基准划分”。

所以可以看到,


成块性 v.s. 有效信息/动力学可逆性

当一个马尔科夫矩阵满足成块性时,同一组的行会受到约束而相似,所以这种成块的矩阵的EI或动力学可逆性也会很高。这说明了成块性和因果涌现之间存在潜在联系。但是,这两个概念并不完全等价。成块性的目标是想在保证宏观和微观动力学的一致性的情况下对状态做分组,而因果涌现的主要目标是想最大化粗粒化之后宏观动力学的因果效应强度。

在普通的成块性例子里,粗粒化之后往往能得到一个 EI 或近似动力学可逆性都很高的宏观动力学。但是,当我们考虑一些极端例子时,如下图所示,我们可以强行让每一组里的状态转移概率尽可能不一样,但加总后又恰好一样,从而仍满足成块性的划分。按照这种成块划分去做粗粒化,并不一定会得到高的 EI 或近似动力学可逆性。相反,用另一个不成块的分组来做粗粒化,得到的宏观动力学 EI 和近似动力学可逆性反而可能更高。这一现象说明:

  1. 基于成块性的最优划分不一定等于基于近似动力学可逆性的最优划分。
  2. 成块性的分组数量是基于马尔科夫链转移矩阵的秩的,而近似动力学可逆性的分组数量是对奇异值谱的截断位置。
成块性和因果涌现.png


划分方法

上述方法基本仅仅给出了优化目标或基本要求,但并不能直接给出粗粒化方案中的划分状态的方法,然而,如何搜索可能的状态空间划分方法以满足或优化上述目标也是一个困难的问题。

当我们想把 [math]\displaystyle{ N }[/math] 个状态的马尔科夫链划分成 [math]\displaystyle{ M }[/math] 组时,待选划分的数量可以由第二类斯特林数(Stirling numbers of the second kind)表示:[math]\displaystyle{ S(N, M) = \frac{1}{M!} \sum_{k=0}^{M} (-1)^{M-k} \binom{M}{k} k^N }[/math]

待选划分的数量极大,比如[math]\displaystyle{ S(20,10)=1.107\times 10^{32} }[/math],几乎不可能通过遍历所有划分来寻找最优解,所以我们需要使用一些方法来有效寻找近似最优解。

逐层贪婪算法

此种方法采取自底向上合并状态(Bottom-Up)的层次聚类(Hierarchical Clustering)方法,即从最细粒度(每个状态单独成组,一共 [math]\displaystyle{ N }[/math] 组)开始,每一步合并两个最相似的组,直到达到目标分组数 [math]\displaystyle{ M }[/math]。

具体步骤如下:

  1. 初始化:每个状态单独成组,初始划分为:[math]\displaystyle{ A_0=\{\{1\},\{2\},…,\{N\}\} }[/math]
  2. 在每一层 [math]\displaystyle{ k }[/math] 内:
    1. 遍历所有两两组对 [math]\displaystyle{ G_p }[/math] 和 [math]\displaystyle{ G_q }[/math],其中[math]\displaystyle{ G_i }[/math]为一状态的集合:
      1. 合并为新组 [math]\displaystyle{ G_{new} = G_p \cup G_q }[/math]
      2. 获得新划分 [math]\displaystyle{ A'_k =A_k \setminus \{ G_p , G_q \} \cup G_{new} }[/math]
      3. 根据此划分获得新动力学(见下文) [math]\displaystyle{ P' }[/math]
      4. 计算 [math]\displaystyle{ P' }[/math] 对应的划分指标(如EI、近似动力学可逆性等)
    2. 更新划分 [math]\displaystyle{ A_{k} }[/math] 为该层最优的划分
  3. 当剩余分组数 [math]\displaystyle{ |A_k| = M }[/math] 时停止

该算法 [19]的复杂度为 [math]\displaystyle{ O \left (\sum_{k=M}^N k^2 \right ) }[/math],介乎[math]\displaystyle{ O \left ( N^2 \right ) }[/math]和[math]\displaystyle{ O \left ( N^3 \right ) }[/math]之间。该算法的缺点在于其容易陷入局部最优。一旦进入了某一分支,就无法跳出。为了找到更好的解,我们可以在其基础上使用遗传算法等方法帮助跳出局部最优。

谱分解

谱分解(Spectral Clustering)的思想,是利用矩阵的谱,也就是特征向量所展示的物理含义来进行聚类和划分。

最大化有效信息的谱分解方法

Griebenow、Klein 和 Erik Hoel[19]先对转移矩阵做特征分解,把状态嵌入到由特征向量张成的谱空间中,再根据谱空间中的几何邻近性来聚类,最后据此构造宏观动力学。它的优点是能把“哪些状态在动力学上相似”转化为一个可操作的聚类问题。

输入:原始包含 [math]\displaystyle{ N }[/math] 个状态的马尔科夫链 [math]\displaystyle{ P }[/math] 和距离超参 [math]\displaystyle{ \epsilon }[/math];

  1. 对 [math]\displaystyle{ P }[/math] 进行矩阵的特征值分解,得到特征值集合 [math]\displaystyle{ \Lambda=\{\lambda_i\}^N_{i=1} }[/math] 与特征向量集合 [math]\displaystyle{ E=\{e_i\}^N_{i=1} }[/math];
  2. 构建有偏特征向量集合 [math]\displaystyle{ E’=\{\lambda_ie_i|\lambda_i\neq 0\} }[/math]。这里用特征值对特征向量加权,并去掉零特征值对应的方向,因为后者更多反映动力学的简并性,而非真正有区分度的大尺度结构;
  3. 依据 [math]\displaystyle{ E' }[/math] 计算状态间的相似矩阵 [math]\displaystyle{ D }[/math]:如果两个状态位于彼此的马尔科夫毯中,就用谱嵌入后的 cosine相似度来衡量其接近程度;否则把它们视为“不可直接聚类”的状态对,赋予极大距离;
  4. 基于相似度矩阵 [math]\displaystyle{ D }[/math] 和超参 [math]\displaystyle{ \epsilon }[/math],使用OPTICS算法在谱空间中聚类,得到候选划分,并从中选出 EI 较高的结果。

忽略零特征值对应的特征向量,实际上是在排除那些不携带有效区分信息的方向;保留下来的谱坐标则更能刻画状态在动力学中的“角色相似性”。该方案的时间复杂度约为 [math]\displaystyle{ O(N^3) }[/math],局限在于特征分解本身成本较高,而且聚类阈值 [math]\displaystyle{ \epsilon }[/math] 往往需要线性搜索或人为调参。

寻找成块划分的谱分解方法

当马尔科夫矩阵存在成块或近成块结构时,谱方法也可以直接用来恢复划分。实际问题中的难点在于:状态顺序常常被打乱,矩阵里还会混入噪声,因此我们往往既不知道它是否严格成块,也不知道应当分成几组。

图7:Zhang[4] 文章中的示意图。图中左面四个矩阵都是成块的马尔科夫矩阵,而右面的P_2是一个噪声矩阵,(P_1)^T P_2 = 0

针对这一问题,Anru Zhang[4]提出先构造稳态流矩阵(frequency matrix)[math]\displaystyle{ F = \mathrm{diag}(\pi) \cdot P }[/math],再利用它的奇异值结构来恢复成块划分。其中 [math]\displaystyle{ \pi }[/math] 是 [math]\displaystyle{ P }[/math] 的稳态分布,[math]\displaystyle{ F_{ij} = \pi_i P_{ij} }[/math] 表示稳态下从状态 [math]\displaystyle{ i }[/math] 流向状态 [math]\displaystyle{ j }[/math] 的长期概率流量。

  1. 获取一个 [math]\displaystyle{ n\times n }[/math] 维的马尔科夫矩阵 [math]\displaystyle{ P }[/math],或者马尔科夫链的频率采样,并由此得到频率矩阵 [math]\displaystyle{ F }[/math];
  2. 获取马尔科夫转移矩阵的秩 [math]\displaystyle{ r }[/math],或者把 [math]\displaystyle{ r }[/math] 视为待选超参;
  3. 对 [math]\displaystyle{ F }[/math] 进行 SVD 分解,[math]\displaystyle{ F = U \Sigma V^T }[/math];
  4. 取前 [math]\displaystyle{ r }[/math] 个右奇异向量,把每个状态映射成向量 [math]\displaystyle{ V_{[i,:r]} }[/math],再在这些点上做聚类。等价地,可以写成寻找下式对应的最优划分 [math]\displaystyle{ \Omega_1,\Omega_2,\ldots,\Omega_r }[/math]:

[math]\displaystyle{ \Omega_1, \Omega_2, ... , \Omega_r = {\arg\min}_{\Omega_1, \Omega_2, ... , \Omega_r} \min_{v_1, v_2, ... , v_r} \sum_{s=1}^r \sum_{i \in \Omega_s} || V_{[i,:r]} - v_s ||_2^2 }[/math]

其中 [math]\displaystyle{ v_s }[/math] 是第 [math]\displaystyle{ s }[/math] 类的中心向量。这个目标函数和 K-Means聚类 本质上是同一类问题:让同一组状态在主奇异子空间中的坐标尽可能靠近。因此,实践中常常直接对这些谱坐标做 K-Means聚类。

为什么是对 [math]\displaystyle{ F }[/math] 做 SVD,而不是直接对 [math]\displaystyle{ P }[/math] 做 SVD?因为 [math]\displaystyle{ P }[/math] 编码的是条件转移分布,而 [math]\displaystyle{ F }[/math] 编码的是稳态下的长期概率流结构。对 [math]\displaystyle{ F }[/math] 做 SVD,相当于在稳态流网络中寻找最主要的流通道,也会自然压低低频状态的干扰。对于成块的情形,同一组状态在主右奇异向量上的坐标会呈现近似分段常数结构,因此能够被较稳定地恢复出来。

这种方法适合用于近成块且存在明显谱间隙的系统:前几个奇异值主要反映组间的大尺度流动,后面的奇异值更多对应组内局部波动。这对应于一种特殊但常见的物理场景:系统各状态群之间的长期概率流量远弱于群内流量,形成类似“亚稳态”的结构,粒子长时间在群内游走,偶尔才跨越群边界。当 [math]\displaystyle{ P }[/math] 不是成块,而是近似成块(near-lumpable)时,可以预期该方法在偏离量较小时仍具有稳定性,但需要满足一个关键条件:组间动力学的奇异值必须主导 [math]\displaystyle{ F }[/math] 的整体奇异值谱。

但需要注意的是,成块性的数学定义本身并不要求组间流量必须强于组内流量。成块性只要求同组状态到各组的转移概率总量相同,并不约束转移的方向强弱。设 [math]\displaystyle{ F = F_1 + F_2 }[/math],其中 [math]\displaystyle{ F_1 }[/math] 是秩为 [math]\displaystyle{ r }[/math] 的成块部分,[math]\displaystyle{ F_2 }[/math] 是组内局部动力学部分,两者行空间正交([math]\displaystyle{ F_1 F_2^T = 0 }[/math]) [math]\displaystyle{ \sigma_r(F_1) \gg \|F_2\| }[/math] 论文给出了分组错误率的上界:

[math]\displaystyle{ \text{Misclassification Rate} \leq \min \left \{ C \frac{\pi_{\max} r \tau_* (\log n)^2/n + \min \left \{ r (\sigma_1(F_2))^2, \|F_2\|_F^2 \right \} }{\sigma_r^2(F_1)}, r \right \} }[/math]

其中 [math]\displaystyle{ C }[/math] 为常数,[math]\displaystyle{ \pi_{max} }[/math] 为稳态分布中最大的的状态概率,[math]\displaystyle{ \tau_* }[/math] 为马尔科夫链的混合时间,[math]\displaystyle{ \sigma_1(F_2) }[/math] 为 [math]\displaystyle{ F_2 }[/math] 的最大奇异值。

分子中第一项是有限样本的统计误差,第二项是近似成块带来的结构偏差,分母是组间动力学的信号强度。因此,当系统虽然满足成块性、但组内转移概率反而很大时(即 [math]\displaystyle{ |F_2| \gtrsim \sigma_r(F_1) }[/math],分母较小而分子较大会拉高错误率上界),SVD 提取的前 [math]\displaystyle{ r }[/math] 个方向会反映组内的局部动力学模式而非组间的成块结构,此时该方法会给出错误的分组。


基于动力学可逆性的谱分解方法

张江等人关于动力学可逆性的文章[18]提供了另一种基于 SVD 的粗粒化方法。它的目标不是恢复某种预设的成块结构,而是把马尔科夫动力学看作一个通信信道,删除马尔科夫动力学中的冗余通路。这里每个状态的出转移分布都可以看作一条信息通路;如果若干行向量彼此线性相关,那么这些通路在信息上就是冗余的,即使它们在概率分布上看起来并不相似。

这也是它和 EI 的一个重要差别:EI 只能度量行向量之间的分布相似性,而动力学可逆性 [math]\displaystyle{ \Gamma_\alpha }[/math] 还能捕捉线性依赖关系。换句话说,两个状态的出流模式可以彼此很不像,但只要在代数上是冗余的,[math]\displaystyle{ \Gamma_\alpha }[/math] 仍能识别出来。算法的基本流程见“近似动力学可逆性”一节。

虽然它与前面的 Zhang & Wang 方法都包含 “SVD + K-Means” 这一步,但两者的目标不同:前者更关心从稳态流结构中恢复成块划分,后者更关心删除 TPM 中的冗余信道。相应地,前者分解的是 [math]\displaystyle{ F=\mathrm{diag}(\pi)P }[/math],因此天然带有稳态权重;后者直接分解 [math]\displaystyle{ P }[/math],更强调转移矩阵本身的线性冗余结构。

其他方法

寻找最优划分本质上是一个组合优化问题,所以我们也能使用其他合适的方法来解决它。

遗传算法

遗传算法是一种受生物进化启发的全局优化算法,通过在每一步迭代过程中模拟自然选择(选择适应度高的个体作为父代)、交叉(通过父代基因重组生成子代)和变异(以一定概率随机改变子代基因)等机制,在解空间中搜索最优解。

梯度下降和深度学习

我们也可以使用梯度下降,搭配图神经网络等深度学习方法来寻找最优划分。比如,在 Klein 和 Hoel 等人的论文[17],以及 Griebenow 等人[19]的文章中,输入原始网络和对应的邻接矩阵。把粗粒化方案写成一个分组矩阵,[math]\displaystyle{ M\in \mathbb{R}^{N×K} }[/math],[math]\displaystyle{ K }[/math]表示分组的大小,其中矩阵里面的每个元素[math]\displaystyle{ m_{iμ}=Pr⁡(v_i\in v_{\mu}) }[/math],并将其参数化,这样就能把最终的粗粒化网络表达为分组矩阵与原始网络邻接矩阵的乘积,[math]\displaystyle{ M^T A M }[/math],并进行归一化,由此计算 EI。在每一步中使用带动量的梯度下降方法优化 [math]\displaystyle{ M }[/math] ,优化目标是最大化宏观网络的有效信息EI。

在上述逐层贪婪算法和谱分解方法的解并未达到最优的情况下,我们可以把它们作为初始解,再借助遗传算法和梯度下降等方法帮助跳出局部最优,从而找到新的全局更优解。


对划分的限制

在上述算法介绍中,并未特别强调各方法对划分节点的限制。但在实际问题场景里,我们往往需要在寻找划分时加入限制,以满足问题本身的特性。

比如,在社区检测问题中,具有只允许相邻的节点聚合在一起的约束。这要求状态划分必须遵循状态转移图的拓扑结构,同一聚合状态内的节点必须在状态转移图中直接相连或通过短路径可达(如K-邻域)。这种限制确保聚合状态具有空间连贯性,同时大幅减少可行划分数量。

在某些因果涌现分组方法[17]中,划分并不是在所有状态之间任意搜索,而是先用马尔科夫毯(Markov blanket)作为结构约束。在概率图中,一个节点的马尔科夫毯通常包括它的父节点、子节点以及子节点的其他父节点。给定马尔科夫毯后,该节点与图中其余节点条件独立;因此,马尔科夫毯可以理解为一个状态与外部系统发生直接依赖关系的局部接口。

在粗粒化问题中,如果两个状态具有相同或高度重合的马尔科夫毯,说明它们面对相似的局部依赖环境,可以作为候选聚合对象。这个约束的作用是缩小划分搜索空间,并提高聚合后宏观状态的局部结构一致性。具有相同马尔科夫毯的状态与其他状态之间的依赖关系相同,可以认为它们在系统中扮演相似的角色。将这些状态聚合成一个宏观状态,相当于识别出了系统中的一个功能模块或模式。

也有许多划分问题不需遵从这两类限制。比如在学术网络、网页流网络的场景上的模式识别时,我们可能会判断某些从未合作过的学者具有相似的行为模式,而把他们划分在一组内。


聚合方法

当我们决定了最优划分之后,就能根据该划分得到聚合后的矩阵。

假设划分是把若干个状态合成一组,而其他状态各自成组,我们要做的是把原矩阵中待聚合状态对应的若干行和若干列,分别聚合成一行和一列,而其他元素保持不动。同时,聚合后的新马尔可夫转移矩阵应满足归一化条件,即所有行和为1。

图8:矩阵聚合

一种最直接的聚合方式是按照下面的两步进行的:把一个 [math]\displaystyle{ N \times N }[/math] 的矩阵聚合成 [math]\displaystyle{ N' \times N' }[/math] 的矩阵,可以分成两步:

  1. 把待聚合的列直接进行相加,得到 [math]\displaystyle{ N \times N' }[/math] 的矩阵;
  2. 把待聚合的行进行加权平均,得到 [math]\displaystyle{ N' \times N' }[/math] 的矩阵。


入流的聚合方法:直接加和

我们可以把马尔科夫链看作一个以状态为节点、以转移概率为边的网络。

我们把待聚合节点都放进一个大节点的框 [math]\displaystyle{ M }[/math] 里,会发现每个外面的节点转移进大节点的概率就等于转移进 [math]\displaystyle{ M }[/math] 中所有节点的概率的和,即[math]\displaystyle{ p_{iM} = \sum_{j \in M} p_{ij} }[/math]。如此就能满足每行加和为一的条件。

图9:入流聚合

所以,对入流(若干列)的聚类方法就是把这若干列直接加和,得到的列等于各节点转移到大节点的概率。


出流的聚合方法:加权平均

由于入流是多对一关系,计算大节点的入流方式比较简单,只有一种:直接加和。而出流是一对多关系,所以计算大节点的出流则相对复杂,需要考虑待聚合节点之间的随机游走和状态分布,来计算流向每个节点的概率。我们可以把出流聚合看作对若干行做线性加权平均。

权重一致的平均

在现有方法中,很多都采用最简单的思路:对状态做权重一致的平均。当我们只看矩阵本身,认为所有状态是平等的,就会给它们分配相同权重,即[math]\displaystyle{ p_{Mi} = \frac{1}{n} \sum_{j \in M} p_{ji} }[/math]。

但是,这种情况只适用于成块的状态划分:由于成块性要求同属一组的微观状态转移到另一组的概率相同,所以这样的平权平均没有问题(成块划分使用任意非负且和为一的权重都会有相同的结果[3][4])。但当我们对不成块的状态划分进行平权平均时,会导致宏观和微观的马尔科夫动力学不一致。

权重为平稳分布的平均

另一种聚合方式是要首先求得马尔可夫链的稳态分布,并将稳态分布下的各连边的稳态概率看作一个流网络,有大量的随机游走粒子在其中流动。当一个随机游走粒子进入大节点时,它下一步游走到大节点内或大节点外其他节点的概率,应当取决于它当前更可能位于大节点内的哪个小节点。

设游走粒子在 [math]\displaystyle{ t }[/math] 时刻处于小节点 [math]\displaystyle{ j }[/math] 的概率为[math]\displaystyle{ s_j(t) }[/math],则聚合节点到外部节点 [math]\displaystyle{ i }[/math] 的下一步概率可写为[math]\displaystyle{ p_{Mi}(t) = \frac{1}{\sum_{j \in M} s_j(t)} \sum_{j \in M} s_j(t)\ p_{ji} }[/math],而留在聚合节点内部的概率为[math]\displaystyle{ p_{MM}(t) = \frac{1}{\sum_{j \in M} s_j(t)} \sum_{j \in M} s_j(t)\sum_{k \in M} p_{jk} }[/math]。

假设游走已经达到平稳分布[math]\displaystyle{ \pi }[/math],游走粒子处于不同小节点 [math]\displaystyle{ j }[/math] 的概率为[math]\displaystyle{ \pi_j }[/math],则聚合节点到外部节点 [math]\displaystyle{ i }[/math] 的概率为[math]\displaystyle{ p_{Mi} = \frac{1}{\sum_{j \in M} \pi_j} \sum_{j \in M} \pi_j\ p_{ji} }[/math],而留在聚合节点内部的概率为[math]\displaystyle{ p_{MM} = \frac{1}{\sum_{j \in M} \pi_j} \sum_{j \in M} \pi_j\sum_{k \in M} p_{jk} }[/math]。按照平稳分布为权重对若干行做加权平均,得到的行仍会满足加和为一的条件。

为了说明如何获得简化的TPM,首先定义一个矩阵,称为稳态流矩阵:[math]\displaystyle{ F_{ij} \equiv \pi_i \cdot P_{ij}, \, \forall i,j \in [1, N], }[/math]。其中,[math]\displaystyle{ \pi }[/math] 是 [math]\displaystyle{ P }[/math] 的稳态分布,即满足[math]\displaystyle{ \pi P = \pi }[/math]。

其次,我们将根据 [math]\displaystyle{ \Phi }[/math] 和 [math]\displaystyle{ F }[/math] 得到粗粒化后的稳态流矩阵:

[math]\displaystyle{ F' = \Phi^T \cdot F \cdot \Phi, }[/math]

其中,[math]\displaystyle{ F' }[/math] 是粗粒化后的平稳流矩阵。在归并的过程中,保持所有矩阵上的流量总量是不变的。

最后,粗粒化后的 [math]\displaystyle{ P' }[/math] 可以通过归一化每一行直接得到:

[math]\displaystyle{ P'_i = \frac{F'_i}{ \sum_{j=1}^{N'} (F'_i)_j }, \, \forall i \in [1, N']. }[/math]

将 [math]\displaystyle{ N' }[/math] 个 [math]\displaystyle{ P'_i }[/math] 拼在一起就能得到最终的 [math]\displaystyle{ P' }[/math]。可以证明,这样的粗粒化方案能够保证粗粒化操作和时间演化操作具有可交换性(证明可以参考[18])。

按照该方法进行加和平均,得到的宏观和微观马尔科夫动力学就会满足一致条件。虽然这个方法直观明了,但它的问题在于:在系统演变到平稳分布之前,需要计算每个 [math]\displaystyle{ t }[/math] 时刻的状态分布概率 [math]\displaystyle{ s_j(t) }[/math]。在达到平稳分布前,直接使用[math]\displaystyle{ \pi }[/math]来聚合仍不会满足一致性,因此操作起来比较麻烦;而 HOM 的计算方式则可以解决这个问题。

高阶网络聚类 Higher-Order Macronodes (HOM)

在 Klein 等人的论文[17],作者采用了高阶依赖项建模(HOMs)[20]来对网络节点进行聚合。其目的是保证聚合后的宏观网络与原始微观网络具有相同的随机游走动力学特性。该方法也可以用于合并马尔科夫链的状态。

Higher-order Macronodes 是由 Jian Xu 等人[20]提出的一种方法:根据给定的微观马尔科夫链和节点划分方案,得出一个与微观转移矩阵保持一致的宏观转移矩阵。下面将简单介绍构建宏观 HOM 转移概率矩阵的计算方法。具体方法请参考网络归并。

假设我们把若干个节点视为待合并的节点[math]\displaystyle{ A }[/math]。

  1. 计算合并节点的入流:对于其他节点到待合并节点[math]\displaystyle{ A }[/math]的连边,直接加和即可。
  2. 计算合并节点的出流:对于待合并节点[math]\displaystyle{ A }[/math]到其他节点的连边,需要考虑三种情况:
    1. 当待合并节点[math]\displaystyle{ A }[/math]间存在连边时;
    2. 当待合并节点[math]\displaystyle{ A }[/math]指向同一个输出节点时;
    3. 当待合并节点[math]\displaystyle{ A }[/math]指向不同输出节点时;
  3. 待合并节点[math]\displaystyle{ A }[/math]的出流的计算方式都不一样,如图(10)所示:
图10:HOM节点边权合并示意图

论文[17]中提到,原始微观动力学与使用该方法聚合后的宏观动力学能够满足其定义的一致性。


本文介绍的粗粒化文章使用上述策略的情况如下:

  1. 任意权重:所有成块性/可交换性文章[3][4][7][8][9]
  2. 平权:Erik Hoel的两篇因果涌现文章[15][16]
  3. 平稳分布权重:Griebenow & Klein & Erik Hoel 在网络上的因果涌现[19]、动力学可逆性[18]、Barnett & Anil Seth的动力学解耦[13]
  4. HOM:Klein & Erik Hoel 另一篇在网络上的因果涌现 [17]

Barnett & Seth 的动力学解耦文章中[13] 把宏观变量定义为确定性粗粒化 [math]\displaystyle{ Y_t=f(X_t) }[/math],实际上用到了平稳分布加权平均。如果我们额外假设微观过程 [math]\displaystyle{ X_t }[/math] 是有限状态一阶 Markov 链,且忽略或不显式保留 [math]\displaystyle{ Y }[/math] 的历史依赖,那么宏观块 [math]\displaystyle{ M }[/math] 的出流应按组内微观状态在平稳分布下的条件概率加权:[math]\displaystyle{ P'_{MN}=\sum_{i\in M} \frac{\pi_i}{\sum_{k\in M}\pi_k} \sum_{j\in N}P_{ij} }[/math]。因此,在这些附加假设下,它可以归为平稳分布加权聚合。

复杂网络中的粗粒化方法

马尔可夫链与复杂网络具有天然的对应关系:一方面,马尔可夫链本身可以看作是一个有向加权的网络,而另一方面,可以针对任何复杂网络定义随机游走动力学,该动力学可以直接对应到一个马尔可夫链。因此,复杂网络的粗粒化方法也可以借鉴作为马尔可夫链的粗粒化方法。

复杂网络中的粗粒化方法非常多,常见思路包括:根据社团结构或模块结构合并节点,根据谱嵌入或扩散距离寻找动力学上相近的节点,根据网络模体或功能单元识别可重复的局部结构,也可以利用复杂网络上的统计物理方法、图神经网络等方法,在保留全局流量、扩散过程或拓扑功能的前提下压缩网络。不同方法关注的对象并不完全相同:有的主要保留拓扑结构,有的保留随机游走或扩散动力学,有的则保留功能模块和多尺度组织。

这些方法与马尔科夫链粗粒化有两层联系。首先,马尔科夫链本身可以看作一种加权有向复杂网络:状态对应节点,转移概率对应有向边权。因此,我们可以先把马尔科夫链的状态转移图视为一个网络,并使用一些有向图方法来生成候选划分。例如,社团检测会倾向于合并转移联系密集、概率流交换频繁的一组状态;谱聚类会根据大尺度连通结构或扩散结构合并状态;角色相似性方法则会把具有相似入流和出流模式的状态放在一起。如果暂时忽略方向、权重或具体转移概率,则可以把马尔科夫链退化为一个普通的无向拓扑网络。这样做会丢掉一部分动力学信息,但也能使用更广泛的复杂网络粗粒化工具,如无向社团检测、模块识别、网络模体识别和层次聚类等。此时得到的划分可作为候选分组,再根据成块性、稳态流、EI、动力学可逆性等指标进行选择。

其次,许多复杂网络上都可以定义随机游走、扩散或传播动力学。一旦规定游走规则,网络就会诱导出一个转移概率矩阵 [math]\displaystyle{ P }[/math],也就是一条马尔科夫链。于是,复杂网络的粗粒化可以转化为马尔科夫链的状态粗粒化:我们不仅关心哪些节点拓扑上相近,也关心哪些节点在随机游走中具有相似的动力学作用。在这个意义上,成块性可以检验合并后的宏观节点是否仍然形成封闭马尔科夫链;稳态流矩阵可以刻画长期概率流通道;谱方法可以寻找扩散尺度上的相似节点;HOM 则可以在节点合并时保留随机游走的高阶路径依赖。因此,马尔科夫链粗粒化工具可以反过来用于复杂网络上的随机游走动力学压缩。

本节只简要介绍几类与马尔科夫链关系较近的复杂网络粗粒化方法,分别对网络上的随机游走动力学、扩散动力学,和多尺度功能模块进行粗粒化。虽然这些方法并不直接给出马尔科夫链的最优粗粒化,但可以为状态划分、网络归并和动力学保真提供有用的启发。

复杂网络的粗粒化方法

复杂网络上的随机游走和谱粗粒化

对一个复杂网络而言,网络本身只描述了节点间的静态关系,并不天然具备某种动力学,所以需要人为定义一套动力学。常见做法是借助随机游走子的概念,在网络上定义一套随机游走动力学。

假设游走子在每一步都会按照边权来进行转移,这种节点间的转移动力学可以用马尔科夫链表示。马尔科夫链的状态对应着复杂网络上的节点,而转移概率对应着网络节点之间的连边权重。

考虑一个连通图 [math]\displaystyle{ G }[/math],邻接矩阵为 [math]\displaystyle{ A }[/math],假设图上有 [math]\displaystyle{ N }[/math] 个节点,分别为 [math]\displaystyle{ 1,2,...,N }[/math],我们定义从节点 [math]\displaystyle{ i\in\left\{1,2,...,N\right\} }[/math] 到节点 [math]\displaystyle{ j\in\left\{1,2,...,N\right\} }[/math] 的转移概率为 [math]\displaystyle{ w_{ij}=\frac{a_{ij}}{\sum_{j}a_{ij}} }[/math](即对边权进行行归一化),其中 [math]\displaystyle{ a_{ij} }[/math] 表示邻接矩阵的第 [math]\displaystyle{ i }[/math] 行第 [math]\displaystyle{ j }[/math] 列元素,对应从节点 [math]\displaystyle{ i }[/math] 到节点 [math]\displaystyle{ j }[/math] 的边的权重。

由于该网络是连通图,所以分母 [math]\displaystyle{ \sum_{j}a_{ij} }[/math] 会大于[math]\displaystyle{ 0 }[/math],而且满足行归一化,即 [math]\displaystyle{ \sum_{j=1}^N {w_{ij}}=1 }[/math]。这样我们就定义了一个马尔科夫链,状态空间为 [math]\displaystyle{ \left\{1,2,...,N\right\} }[/math],从状态 [math]\displaystyle{ i }[/math] 到状态 [math]\displaystyle{ j }[/math] 的转移概率为 [math]\displaystyle{ w_{ij} }[/math],转移矩阵为 [math]\displaystyle{ W=\left\{w_{ij}\right\}_{N \times N} }[/math](详情请见复杂网络中的因果涌现词条)。

为了对网络的随机游走矩阵 [math]\displaystyle{ W_{ij} = \frac{A_{ij}}{\sum_k A_{kj}} }[/math] 进行降维与简化,研究者提出了一种基于谱理论的粗粒化方法[21]。这里采用的是列归一化写法:每一列的和为 1,因此指标 [math]\displaystyle{ j }[/math] 表示出发节点,[math]\displaystyle{ i }[/math] 表示到达节点。该方法的核心在于选取随机矩阵前 [math]\displaystyle{ S }[/math] 个最大非平凡特征值 [math]\displaystyle{ \lambda^2, \lambda^3, \ldots, \lambda^{S+1} }[/math] 对应的左特征向量 [math]\displaystyle{ \langle u^2 |, \langle u^3 |, \ldots, \langle u^{S+1} | }[/math],将网络节点映射到这些特征向量张成的特征空间。在此基础上,将每个特征向量的值域划分为 [math]\displaystyle{ I }[/math] 个区间,在所有 [math]\displaystyle{ S }[/math] 个特征向量上都落入相同区间的节点会被合并为一个超级节点。

从概率角度看,这一粗粒化过程可以理解为:新超级节点的出边概率由原节点的出边概率加权平均得到,其流入概率则由原节点的流入概率直接相加,其他未合并节点的转移概率保持不变。


结合复杂网络配分函数与图机器学习的网络流方法

本节介绍的方法主要参考文章[22],该方法结合了复杂网络上的统计物理方法与图神经网络方法,从而对复杂网络进行粗粒化。其核心思想是,粗粒化前后两个网络的配分函数应该保持不变,将二者的误差作为损失函数,训练一个图神经网络,从而提炼出如何将大网络转化为粗粒化的小网络的映射方法。之所以将复杂网络的配分函数作为一个重要的评价指标,是因为该函数不仅可以体现网络的拓扑结构特征,也能够反映网络在多个时间尺度上的动力学特征,因而反映整体网络功能。实验证明,这套粗粒化网络的方法可以应用到大量的真实网络中,并得到了良好的效果[22]。

下面介绍具体的方法。首先,我们用扩散方程 [math]\displaystyle{ \frac{d\psi}{d\tau}=-L\psi }[/math] 来描述网络的扩散过程,其中拉普拉斯矩阵 [math]\displaystyle{ L=D-A }[/math] 描述网络的拓扑结构,[math]\displaystyle{ A }[/math] 是邻接矩阵,[math]\displaystyle{ D }[/math] 是度矩阵。扩散方程的解为 [math]\displaystyle{ \psi_{\tau}=e^{-\tau L}\psi_0 }[/math],体现了信息随时间的演化,其中 [math]\displaystyle{ e^{-\tau L} }[/math] 为时间演化算子,[math]\displaystyle{ [e^{-\tau L}]_{ij} }[/math] 表示经过时间 [math]\displaystyle{ \tau }[/math] 后从节点 [math]\displaystyle{ j }[/math] 到节点 [math]\displaystyle{ i }[/math] 的流量。

其次,利用复杂网络的配分函数 [math]\displaystyle{ Z_{\tau}=Tr(e^{-\tau L}) }[/math],我们能量化矩阵的全局信息量[22]。矩阵 [math]\displaystyle{ e^{-\tau L} }[/math] 的每个对角元 [math]\displaystyle{ [e^{-\tau L}]_{ii} }[/math] 是第 [math]\displaystyle{ i }[/math] 个节点的信息保留量,对所有对角元求和则能得到所有节点的信息保留量。

通过图神经网络(GNN)的消息传递机制,得到每个节点的向量表征,对节点向量应用Softmax函数,得到分组矩阵 [math]\displaystyle{ S_{N\times M} }[/math],其中元素 [math]\displaystyle{ s_{ij} }[/math] 表示节点 [math]\displaystyle{ i }[/math] 分到第 [math]\displaystyle{ j }[/math] 个超级节点的概率。

通过 [math]\displaystyle{ \hat{A}=\frac{M}{N}S^{T}AS }[/math] 得到宏观邻接矩阵 [math]\displaystyle{ \hat{A}_{M\times M} }[/math],其中 [math]\displaystyle{ A_{N\times N} }[/math] 原始网络的邻接矩阵。归一化配分函数 [math]\displaystyle{ \frac{Z}{N} }[/math] 的平均绝对误差(MAE)作为训练图神经网络的损失函数,通过梯度下降优化GNN的参数。

网络规模从 [math]\displaystyle{ N }[/math] 压缩到 [math]\displaystyle{ N^{'} }[/math] 后,得到新的拉普拉斯矩阵 [math]\displaystyle{ L^{'} }[/math] 和配分函数 [math]\displaystyle{ Z^{'} }[/math],原始网络和压缩后的网络的配分函数的关系是 [math]\displaystyle{ Z=\gamma Z^{'} }[/math],并且压缩后的配分函数和原始网络的配分函数要尽可能相似,不随时间变化而且误差最小。当 [math]\displaystyle{ \tau=0 }[/math] 时,[math]\displaystyle{ Z=N,Z^{'}=N^{'} }[/math],为了使之成立,[math]\displaystyle{ \gamma=\frac{N}{N^{'}} }[/math],所以乘数的有效值为 [math]\displaystyle{ \frac{N}{N^{'}} }[/math]。

马尔科夫链是描述状态仅依赖当前状态的随机过程,通过拉普拉斯矩阵与转移概率的关系,将上述方法应用到马尔科夫链中,在离散时间马尔科夫链中,转移概率矩阵为 [math]\displaystyle{ P }[/math];在连续时间马尔科夫链中,[math]\displaystyle{ P(t)=e^{-tL} }[/math] 是经过时间 [math]\displaystyle{ t }[/math] 的转移概率矩阵。配分函数 [math]\displaystyle{ Z_t=Tr(e^{-tL}) }[/math] 可以理解为马尔科夫链的所有状态的信息保留总和。

进行稳态分析的时候,当 [math]\displaystyle{ t\rightarrow \infty }[/math] 时,[math]\displaystyle{ e^{-tL} }[/math] 的极限与稳态分布有关,配分函数的渐近行为反映链的长期信息保留特性。通过分析 [math]\displaystyle{ e^{-tL} }[/math] 的特征向量来识别马尔科夫链的状态聚类。

具有自非相似性的复杂网络的粗粒化方法

文献[23]介绍的方法将复杂网络抽象出来一个“自非相似性”的特性,从而对复杂网络进行粗粒化。自非相似性是指复杂网络的不同粗粒化层面,粗粒化单元和网络模体与整体网络的拓扑结构存在显著差异,而不是保持部分与整体的相似性。不同粗粒化层级的粗粒化单元和网络模体通过功能分化实现整体功能。只要网络存在功能模块的重复级联,这个方法就可以应用,目前该方法在电子电路(ISCAS89基准电路)和生物信号网络(哺乳动物蛋白质相互作用网络)成功应用,。

定义粗粒化单元(CGU):它是最适合作为网络粗粒化节点的模式,具有明确输入、输出端口和内部计算节点,可视为一种基本电路组件。CGU 的集合构成“元素字典”,粗粒化网络则可视为由该字典构建的版本。部分节点被 CGU 替换,从而实现网络压缩。

CGU的评价维度(优化的属性):简洁性、简单性和覆盖度。

  1. 简洁性:字典中 CGU 的总数尽可能少。
  2. 简单性:单个 CGU 的结构尽可能简单,由黑箱数 [math]\displaystyle{ H }[/math] 决定,[math]\displaystyle{ H }[/math] 越少越简单。 [math]\displaystyle{ H=I+O+2M }[/math], [math]\displaystyle{ I }[/math] 为输入节点,[math]\displaystyle{ O }[/math] 为输出节点,[math]\displaystyle{ M }[/math] 为混合节点。
  3. 覆盖度:粗粒化后消除的节点和边数尽可能多,被CGU覆盖的节点数 [math]\displaystyle{ P_{covered} }[/math] 和边数 [math]\displaystyle{ E_{covered} }[/math] 尽可能多。

针对CGU的评分函数:[math]\displaystyle{ S=[E_{covered}+\alpha P_{covered}]-[\alpha \sum_{i=1}^{N}n_i H_i+\beta N+\gamma\sum_{i=1}^{N} T_i] }[/math]

  1. 正向项:对应覆盖度,覆盖度高,得分越高。
  2. 负向项:对应简洁性和简单性,CGU 类型越少,内部节点越少,得分越高。

最大化 S 倾向于选择“高频出现、内部节点多且混合节点少” 的 CGU。

采用模拟退火法寻找粗粒化最优的 CGU 集合,仅以网络模体(在原始网络中出现频率显著高于随机网络的子图)作为候选 CGU。 候选CGU的获取方式:先检测 3-6 节点的所有网络模体,再根据每个模体出现的连接特征向量,同一模体不同特征向量的出现视为不同 CGU 候选。

找到最优的粗粒化CGU集合,通过将每个CGU的每一次出现替换为一个节点,形成粗粒化网络。算法可在粗粒化网络上重复执行,获得更高层级的粗粒化。


应用

马尔科夫链粗粒化的应用工作虽然来自不同领域,但结构和思路都非常相似:

  1. 先有一个微观马尔科夫链;
  2. 微观状态数量很多,或者观测上太细;
  3. 根据动力学相似性、观测等价性、研究目标或系统结构,把微观状态合并;
  4. 得到更少的宏观状态;
  5. 再分析粗粒化是否保留了原始系统的重要信息。

需要注意的是,粗粒化的结果不一定仍然是严格的马尔科夫链。很多应用正是通过粗粒化来研究:哪些状态可以合并?合并以后会不会丢失马尔科夫性?如果丢失,能不能构造一个近似的宏观马尔科夫模型?


分子构象动力学

代表文献: Kube & Weber, 2007, A coarse graining method for the identification of transition rates between molecular conformations. [24]

在分子动力学中,蛋白质、小分子或药物分子会不断改变自己的空间形状。例如,蛋白质可能会折叠、展开,或者在几种相对稳定的形状之间切换。这里的“构象”指的是分子的空间形状:同一个分子中,原子连接关系没有变,但整体姿态、弯曲方式或折叠方式不同,就对应不同构象。

分子动力学模拟会产生大量构象。为了建立马尔科夫链,研究者通常先把连续的构象空间离散化:把相近的分子形状归成许多小的构象簇,每个构象簇看作一个微观状态。然后,根据模拟轨迹估计转移概率:如果分子现在处在微观构象 i,过一段时间后转到微观构象 j 的概率是多少?这些概率组成转移矩阵 P,于是分子在不同构象簇之间的随机跳转就可以近似看作马尔科夫链。

粗粒化的目的,是把大量细小构象合并成少数几个有物理意义的“大状态”。例如,一个蛋白质可能有展开态、半折叠中间态和完全折叠态。每个大状态内部仍包含许多细微不同的构象。这类粗粒化主要依据转移动力学,而不是单纯依据几何形状。如果一批微观构象之间经常互相转移,但很少跳到外面,那么它们可以合并成一个宏观状态,称为“亚稳态构象”。“亚稳态”可以理解为:系统会在这个状态附近停留较久,但最终仍可能跳到别的状态。Kube 和 Weber 的工作使用PCCA+方法。该方法通过分析转移矩阵,找出在较长时间内不容易离开的状态集合。粗粒化之后,研究者可以估计这些大状态之间的转移速率,例如蛋白质从展开态变成折叠态大约需要多久。

生化反应网络

合并复杂的分子组合

代表文献:Ganguly, Petrov & Koeppl, 2014, Markov chain aggregation and its applications to combinatorial reaction networks[25];Petrov, 2019, Markov chain aggregation and its application to rule-based modeling.[26]

细胞内部有大量化学反应。一个蛋白质可以被磷酸化,可以和另一个蛋白质结合,也可以从复合物中解离。由于分子可能有多个结合位点和修饰位点,所有可能的分子组合会非常多,即使反应规则本身不复杂,状态数也可能爆炸式增长。在这类模型中,微观状态是系统在某一时刻的完整化学状态。例如,一个微观状态需要说明:系统里有哪些分子、每种分子有多少个、哪些分子彼此结合、哪些位点被修饰等。每一次反应都会把系统从一个微观状态带到另一个微观状态。如果每个反应都有发生速率,那么整个系统就形成一个连续时间马尔科夫链。

完整微观链通常太大,难以直接计算。例如,要求未来某个时间系统处于某类状态的概率,或者求长期稳态分布,都可能因为状态数太多而不可行。同时,研究者往往并不关心所有细节,而只关心某些关键变量,例如某个蛋白是否被激活。因此,可以把对目标变量没有影响或近似没有影响的细节合并掉。粗粒化依据可以是研究目标、分子结构、反应网络对称性,或者成块性条件。粗粒化后得到的宏观状态可能是:“关键蛋白未激活”、“关键蛋白已激活但未结合配体”、“关键蛋白已激活且结合配体”等。

Ganguly 等人的工作研究的是:给定一个连续时间马尔科夫链和一个状态分组,在什么条件下可以在分组后的状态上定义合理的聚合模型。Petrov 的工作进一步讨论如何在 rule-based biochemical models 中自动寻找这种聚合结构。这里的 rule-based model 指的是用反应规则描述系统,而不是一开始就枚举所有可能的微观状态,因此特别适合处理状态空间极大的生化系统。

合并物种集合

代表文献:Mu et al., 2018, Organisation-Oriented Coarse Graining and Refinement of Stochastic Reaction Networks.[27]

另外,在一些化学反应系统中,研究者关心的不是每一种分子此刻具体有多少个,而是系统最终会维持哪些物种。这里的“物种”可以理解为反应网络中的化学成分,例如某种分子、离子或复合物。“chemical organisation” 可以翻译为“化学组织”,指一组在反应中保持封闭、并且有可能自我维持的物种集合。通俗地说,如果系统进入某个化学组织,其中的物种可以通过反应互相支持,不会立刻消失,也不会不断产生完全不属于这个集合的外部物种。

这类模型的微观状态通常是每种物种的分子数。例如,若系统中有三种物种 A,B,C,一个微观状态可以写成[math]\displaystyle{ (n_A,n_B,n_C) }[/math]。每次反应都会改变这些分子数,例如[math]A+B\rightarrow C[/math] 会减少 A 和 B,增加 C。若反应随机发生且各自有发生速率,这些分子数状态之间的跳转就形成一个连续时间马尔科夫链。粗粒化的依据不是“分子数是否接近”,而是当前状态中哪些物种存在,以及这些物种是否形成某个 chemical organisation。粗粒化后得到的宏观状态是不同的化学组织,例如:

  • 只有物种 (A) 能维持;
  • 物种 (A) 和 (B) 能共同维持;
  • 物种 (A,B,C) 形成更大的自我维持结构。

这类粗粒化可以帮助研究者理解反应网络的长期行为,而不必追踪每一个具体分子数状态。

网络社团检测

代表文献:Piccardi, 2011, Finding and Testing Network Communities by Lumped Markov Chains.[28]

网络社团检测的目标是找出网络中联系紧密的节点集合。例如,在社交网络中,社团可能是一群经常互动的人;在网页网络中,社团可能是一组主题相近的网站;在交通网络中,社团可能是一组流动联系紧密的区域。

在这个应用中,微观状态是网络节点。给定一个网络,可以定义一个随机游走过程:如果当前在节点 i,下一步就按照边的权重随机移动到某个邻居节点 j。这样就得到一个节点层面的马尔科夫链。粗粒化的目标是把节点合并成较少的社团。其依据是随机游走的停留性:如果一个节点集合内部连接很强,随机游走进入后不容易离开,那么这个节点集合就可以看作一个社团。这个例子说明,马尔科夫链粗粒化可以为“什么是社团”提供一个动力学解释:社团就是随机游走不容易离开的地方。


城市交通与出租车流

代表文献:Zhang & Wang, 2020, Spectral State Compression of Markov Processes [4];Duan, Ke & Wang, 2019, *State Aggregation Learning from Markov Transition Data.[6]

城市中的出租车、共享单车、手机信令、通勤和迁徙数据都可以看成地点之间的移动记录。例如,一次出租车行程从地点 i 到地点 j,就表示一次从 i 到 j 的转移。如果把城市划分成网格、街区或交通小区,就可以用马尔科夫链描述人在城市空间中的移动。

在这个应用中,微观状态是城市中的细空间单元,例如网格、街区、交通小区、出租车上车/下车区域等。根据大量出行记录,可以估计转移概率,从而得到城市移动马尔科夫链。城市空间单元很多,原始转移矩阵很大。城市空间单元很多,原始转移矩阵很大,直接查看每个网格到每个网格的流动通常很难解释。粗粒化可以帮助识别更大的城市功能区域或交通模式,例如商业区、居住区、旅游区、交通枢纽等。

粗粒化依据是区域之间的移动模式是否相似。如果两个空间单元的人流来源和去向相似,它们可能属于同一类区域。Duan 等人的方法进一步允许“软聚合”:一个区域不一定只属于一个宏观类,而是可以以不同概率属于多个 meta-states。这里的 meta-state 可以理解为从数据中学习出来的宏观状态或潜在类型。这很符合城市空间的实际情况。例如,一个区域白天可能像办公区,晚上可能像餐饮娱乐区,因此它未必只属于一种类型。

因此,粗粒化后得到的是空间单元的聚合。硬分组时,每个网格或街区属于一个宏观区域。软分组时,每个区域对多个宏观类型都有一定隶属度。例如,在曼哈顿出租车数据中,粗粒化后的状态可能对应时代广场一类的高活动商业娱乐区、博物馆或公园附近的旅游区域、World Trade Center 附近的商务区域、交通换乘密集区域。


海洋流与空间流动

代表文献:Faccin, Schaub & Delvenne, 2021, State aggregations in Markov chains and block models of networks. [29]

在海洋、空气污染、生态扩散等系统中,研究者常关心物质如何在空间中移动。例如,海水、漂浮物、污染物或生物幼体可能随着流场移动。如果把空间切成许多小格子,就可以描述一个粒子从一个格子移动到另一个格子的概率。

在这个应用中,微观状态是空间小格子或网络节点。例如,在海洋表层流中,可以把海面划分成许多网格。一个漂浮粒子现在位于网格 i,过一段时间后移动到网格 j 的概率,就是马尔科夫链中的转移概率。细网格数量很多,而研究者通常关心更大尺度的流动结构。例如,哪些区域之间交换频繁?哪些区域在长时间内相对封闭?哪些区域共同形成一个大尺度循环?粗粒化可以把大量空间网格合并成少数动力学模块。

Faccin 等人的方法依据的是“保留长期转移信息”。简单说,一个好的分组应该让我们即使只知道粒子所在的宏观区域,也仍然能尽量预测它在较长时间后会到哪里。该方法用信息论中的互信息来衡量这种保留程度。而粗粒化后得到的分组是空间上的动力学模块。在海洋流中,这些模块可能对应一些大尺度相干区域。粒子在同一模块内部更容易一起移动或长期停留,而不同模块之间的交换相对较少。


Agent-based model 与观点动力学

代表文献:Banisch, Lima & Araújo, 2012, Agent Based Models and Opinion Dynamics as Markov Chains.[30]

Agent-based model是一种个体层面的建模方法。模型中有许多个体,每个个体按照简单规则行动或改变状态。例如,在观点动力学中,每个人可能支持观点 0 或观点 1,并受到邻居影响而改变观点。虽然个体规则简单,但所有个体组合起来以后,系统状态会非常多。

在这类模型中,微观状态是所有个体状态的完整配置。比如,有 N 个人,每个人都有两种观点 0 或 1,那么整个系统的一个微观状态就是一个长度为 N 的 0/1 序列,例如 (1,0,0,1,1,...)。如果每一步随机选择一个人更新观点,系统就会从一个配置随机跳到另一个配置。因此,所有个体配置构成一个马尔科夫链。

研究者通常不关心每个人的具体观点排列,而是关心宏观变量,例如支持观点 1 的人数、两种观点是否达成共识、是否出现极化等。粗粒化的依据就是这些宏观变量。例如,可以把所有“有 k 个人支持观点 1”的微观配置合并在一起。如果系统的更新规则只依赖于这个人数,那么粗粒化后的过程可能仍然是马尔科夫链。如果更新规则还依赖于谁和谁相连、支持者分布在哪些位置,那么只知道人数可能不够,宏观过程就可能不再是马尔科夫链。这个例子说明,马尔科夫链粗粒化可以帮助我们从个体规则推导宏观社会变量的演化;但当宏观变量选得太粗时,也可能丢失关键结构信息。

排队网络与计算机系统

代表文献:Song & Takahashi, 1991, Aggregate Approximation for Tandem Queueing Systems with Production Blocking[31] ;Gambin, Krzyżanowski & Pokarowski, 2008, Aggregation algorithms for perturbed Markov chains with applications to networks modeling[32]

排队系统和计算机系统经常可以用马尔科夫链建模。例如,一个工厂生产线有多个工位,每个工位前面可能有队列;一个服务器系统有多个任务等待处理。系统的完整状态需要记录每个队列有多少任务、每个服务器是否忙碌、是否发生阻塞等。

在这类模型中,微观状态是系统的完整配置。例如,在一个有三个队列的系统中,状态 (2,5,1) 表示三个队列中分别有 2、5、1 个任务。服务完成、任务到达或阻塞解除都会使系统从一个状态跳到另一个状态。因此,排队系统的随机演化可以表示为离散时间或连续时间马尔科夫链。

完整状态空间往往很大。队列越多、容量越大,可能状态数就越多。直接求稳态分布、平均等待时间或阻塞概率会变得很困难。工程上通常更关心平均队长、平均等待时间、系统吞吐量、阻塞概率、服务器利用率等聚合指标。因此,可以通过粗粒化降低模型规模。这类粗粒化通常依据系统结构和耦合强弱。如果一组微观状态之间转移频繁,而与其他组之间联系较弱,就可以把它们看成一个聚合状态。粗粒化后的分组可以理解为聚合性能状态,例如“低拥堵状态”“中等拥堵状态”或“高拥堵状态”。这类粗粒化的主要目的通常不是发现新的自然类别,而是让大规模马尔科夫模型可以近似计算。


从上述工作中可见,在不同领域中,粗粒化的含义可以不同:

  1. 在分子动力学中,粗粒化得到的是几个稳定的分子形状;
  2. 在生化反应网络中,粗粒化用于避免状态空间爆炸;
  3. 在化学反应网络长期行为分析中,粗粒化得到的是能够自我维持的物种集合;
  4. 在网络科学中,粗粒化得到社团;
  5. 在交通和迁徙数据中,粗粒化得到功能区、交通模式或城市类型;
  6. 在海洋流和空间流动中,粗粒化得到大尺度动力学模块;
  7. 在观点动力学中,粗粒化得到宏观社会变量;
  8. 在排队系统中,粗粒化得到可计算的性能状态。
应用领域 代表问题 微观马尔科夫链 宏观状态或近似宏观链 粗粒化依据
分子构象动力学 蛋白质、小分子在不同空间形状之间切换 微观状态是离散化后的分子构象簇;转移矩阵由分子动力学或蒙特卡罗轨迹估计得到 亚稳态构象,例如展开态、中间态、折叠态;可进一步构造构象之间的近似宏观转移矩阵 微观构象之间的转移动力学相似性,尤其是是否形成长期停留的亚稳态集合
生化反应网络 复杂分子组合、蛋白结合和修饰导致状态空间爆炸 微观状态是完整化学状态,包括分子数量、结合关系、修饰状态等;反应速率定义连续时间马尔科夫链 聚合后的化学状态,例如关键蛋白未激活或已激活等;在满足条件时可得到聚合 CTMC 研究目标、反应网络结构、分子对称性或成块性条件
化学反应网络长期行为 分析随机反应网络最终能维持哪些物种 微观状态是各物种的分子数向量,例如 [math]\displaystyle{ (n_A,n_B,n_C) }[/math];随机反应定义连续时间马尔科夫链 化学组织,即能够封闭并自我维持的物种集合;宏观过程描述系统在不同组织之间的转移或趋向 当前状态中哪些物种存在,以及这些物种是否形成化学组织
网络社团检测 在复杂网络中识别联系紧密的节点集合 微观状态是网络节点;随机游走定义节点层面的离散时间马尔科夫链 网络社团;宏观过程可理解为社团之间的随机游走 随机游走在节点集合内的停留性,即进入后是否不容易离开
城市交通与出租车流 从出租车、通勤、迁徙等数据中识别城市功能区或流动模式 微观状态是网格、街区、交通小区等细空间单元;出行记录估计空间单元之间的转移矩阵 功能区、交通模式或城市类型;宏观过程描述这些区域或类型之间的流动 空间单元的人流来源和去向是否相似;也可使用软聚合学习潜在 meta-states
海洋流与空间流动 分析海水、污染物、漂浮物等在空间中的长期输运结构 微观状态是空间网格或网络节点;粒子从一个网格移动到另一个网格的概率定义马尔科夫链 大尺度动力学模块或相干区域;宏观过程描述模块之间的长期输运关系 是否保留长期转移信息,例如相隔多步后的互信息
Agent-based model 与观点动力学 从个体互动规则推导宏观社会状态变化 微观状态是所有个体状态的完整配置,例如长度为 [math]\displaystyle{ N }[/math] 的 0/1 观点序列;随机更新规则定义马尔科夫链 宏观社会状态,例如“有 [math]\displaystyle{ k }[/math] 个人支持观点 1”;是否仍为马尔科夫过程取决于宏观变量是否保留足够信息 选定的宏观变量,例如支持某观点的人数、群体比例、是否达成共识
排队网络与计算机系统 近似计算大型排队系统、生产系统或计算机系统的性能指标 微观状态是完整系统配置,例如各队列长度、服务器状态、阻塞状态等;事件发生定义 DTMC 或 CTMC 聚合性能状态,例如低拥堵、中等拥堵、高拥堵状态;宏观模型用于近似稳态分布和性能指标 系统结构、耦合强弱、状态之间的快速内部转移和弱外部转移


因此,马尔科夫链粗粒化并不是简单地“减少状态数”,而是在减少状态数的同时,尽量保留对研究问题重要的转移动力学。不同应用中的“重要信息”并不相同:它可以是亚稳态结构、目标变量、社团保持性、长期流动信息、宏观社会变量,也可以是系统性能指标。

参考文献

  1. ↑ Hafner, Danijar, et al. "Mastering atari with discrete world models." arXiv preprint arXiv:2010.02193 (2020).
  2. ↑ 2.0 2.1 2.2 Coarse graining. Encyclopedia of Mathematics. URL: http://encyclopediaofmath.org/index.php?title=Coarse_graining&oldid=16170
  3. ↑ 3.0 3.1 3.2 3.3 3.4 Buchholz, Peter. "Exact and ordinary lumpability in finite Markov chains." Journal of applied probability 31.1 (1994): 59-75.
  4. ↑ 4.0 4.1 4.2 4.3 4.4 4.5 Zhang, Anru, and Mengdi Wang. "Spectral state compression of markov processes." IEEE transactions on information theory 66.5 (2019): 3202-3231.
  5. ↑ Singh, Satinder, Tommi Jaakkola, and Michael Jordan. "Reinforcement learning with soft state aggregation." Advances in neural information processing systems 7 (1994).
  6. ↑ 6.0 6.1 Duan, Yaqi, Tracy Ke, and Mengdi Wang. "State aggregation learning from markov transition data." Advances in Neural Information Processing Systems 32 (2019).
  7. ↑ 7.0 7.1 7.2 7.3 7.4 Kemeny, John G., and J. Laurie Snell. Finite markov chains. Vol. 26. Princeton, NJ: van Nostrand, 1969. https://www.math.pku.edu.cn/teachers/yaoy/Fall2011/Kemeny-Snell_Chapter6.3-4.pdf
  8. ↑ 8.0 8.1 Barreto, André MS, and Marcelo D. Fragoso. "Lumping the states of a finite Markov chain through stochastic factorization." IFAC Proceedings Volumes 44.1 (2011): 4206-4211.
  9. ↑ 9.0 9.1 Otsuka, Jun, and Hayato Saigo. "Process theory of causality: a category-theoretic perspective." Behaviormetrika 51.1 (2024): 21-36.
  10. ↑ Zennaro, F. M., Turrini, P., & Damoulas, T. (2023). Quantifying consistency and information loss for causal abstraction learning. arXiv preprint arXiv:2305.04357.
  11. ↑ Beckers, S., Eberhardt, F., & Halpern, J. Y. (2020, August). Approximate causal abstractions. In Uncertainty in artificial intelligence (pp. 606-615). PMLR.
  12. ↑ 12.0 12.1 12.2 Rosas, F. E., Geiger, B. C., Luppi, A. I., Seth, A. K., Polani, D., Gastpar, M., & Mediano, P. A. (2024). Software in the natural world: A computational approach to hierarchical emergence. arXiv preprint arXiv:2402.09090.
  13. ↑ 13.0 13.1 13.2 Barnett, Lionel, and Anil K. Seth. "Dynamical independence: discovering emergent macroscopic processes in complex dynamical systems." Physical Review E 108.1 (2023): 014304.
  14. ↑ B. C. Geiger and C. Temmel, “Lumpings of Markov chains, entropy rate preservation, and higher-order lumpability,” J. Appl. Probab., vol. 51, no. 4, pp. 1114–1132, Dec. 2014, extended version: arXiv:1212.4375.
  15. ↑ 15.0 15.1 15.2 Hoel, Erik P.; Albantakis, L.; Tononi, G. (2013). "Quantifying causal emergence shows that macro can beat micro". Proceedings of the National Academy of Sciences. 110 (49): 19790–19795.
  16. ↑ 16.0 16.1 Hoel E P. When the map is better than the territory[J]. Entropy, 2017, 19(5): 188.
  17. ↑ 17.0 17.1 17.2 17.3 17.4 17.5 Klein B, Hoel E. The emergence of informative higher scales in complex networks[J]. Complexity, 2020, 20201-12.
  18. ↑ 18.0 18.1 18.2 18.3 18.4 18.5 Zhang, J., Tao, R., Leong, K.H. et al. Dynamical reversibility and a new theory of causal emergence based on SVD. npj Complex 2, 3 (2025).
  19. ↑ 19.0 19.1 19.2 19.3 Griebenow R, Klein B, Hoel E. Finding the right scale of a network: efficient identification of causal emergence through spectral clustering[J]. arXiv preprint arXiv:190807565, 2019.
  20. ↑ 20.0 20.1 Xu, J., Wickramarathne, T. L., & Chawla, N. V. Representing higher-order dependencies in networks[J]. Science advances, 2016, 2(5), e1600028.
  21. ↑ Gfeller D, De Los Rios P. Spectral coarse graining of complex networks[J]. Physical review letters, 2007, 99(3): 038701.
  22. ↑ 22.0 22.1 22.2 Zhang, Zhang, et al. "Coarse-graining network flow through statistical physics and machine learning." Nature Communications 16.1 (2025): 1605.
  23. ↑ Itzkovitz, Shalev, et al. "Coarse-graining and self-dissimilarity of complex networks." Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 71.1 (2005): 016127.
  24. ↑ Kube, Susanna, and Marcus Weber. "A coarse graining method for the identification of transition rates between molecular conformations." The Journal of chemical physics 126.2 (2007).
  25. ↑ Ganguly, Arnab, Tatjana Petrov, and Heinz Koeppl. "Markov chain aggregation and its applications to combinatorial reaction networks." Journal of mathematical biology 69.3 (2014): 767-797.
  26. ↑ Petrov, Tatjana. "Markov chain aggregation and its application to rule-based modelling." Modeling Biomolecular Site Dynamics: Methods and Protocols. New York, NY: Springer New York, 2019. 297-313.
  27. ↑ Mu, Chunyan, et al. "Organisation-oriented coarse graining and refinement of stochastic reaction networks." IEEE/ACM Transactions on Computational Biology and Bioinformatics 15.4 (2018): 1152-1166.
  28. ↑ Piccardi, Carlo. "Finding and testing network communities by lumped Markov chains." PloS one 6.11 (2011): e27028.
  29. ↑ Faccin, Mauro, Michael T. Schaub, and Jean-Charles Delvenne. "State aggregations in Markov chains and block models of networks." Physical Review Letters 127.7 (2021): 078301.
  30. ↑ Banisch, Sven, Ricardo Lima, and Tanya Araújo. "Agent based models and opinion dynamics as Markov chains." Social Networks 34.4 (2012): 549-561.
  31. ↑ Song, Yu, and Yukio Takahashi. "Aggregate approximation for tandem queueing systems with production blocking." Journal of the Operations Research Society of Japan 34.3 (1991): 329-353.
  32. ↑ Gambin, Anna, Piotr Krzyżanowski, and Piotr Pokarowski. "Aggregation algorithms for perturbed Markov chains with applications to networks modeling." SIAM Journal on Scientific Computing 31.1 (2008): 45-73.