马尔科夫链的粗粒化

马尔科夫矩阵是指满足每一行和为1的条件的方阵，而马尔科夫链指的是一个n维的状态的序列[math]\displaystyle{ s^{(t)}\ = \{1, ..., n\}^{(t)} }[/math]，每一步的状态转换都有马尔科夫矩阵[math]\displaystyle{ P }[/math]决定，即[math]\displaystyle{ p(s^{(t+1)}) = p(s^{(t)}) P }[/math]。

[math]\displaystyle{ P }[/math]的每一行对应的每个状态转移到其他状态的概率，即[math]\displaystyle{ p_{ij} = p(s^{(t+1)} = j | s^{(t)} = i) }[/math]。比如，当[math]\displaystyle{ s^{(t)} }[/math]等于第一个状态的时候，[math]\displaystyle{ P }[/math]的第一行展示了[math]\displaystyle{ s^{(t+1)} }[/math]状态的概率。

对马尔科夫链做粗粒化的意义是什么呢？我们看到文献中着重强调这两点：

1. 我们在研究一个超大的系统的时候，比如复杂城市系统时，并不会关注每一个微观的状态的变化，在粗粒化中我们希望能过滤掉一些我们不感兴趣的噪声和异质性，而从中微观尺度中总结一些中尺度或宏观尺度的规律；
2. 有些状态的转移概率非常相似，所以可以被看成同一类状态，对这种马尔科夫链做partitioning可以减少系统表示的冗余性；
3. 在用到马尔科夫决策过程的强化学习里，对马尔科夫链做粗粒化可以减少状态空间的大小，提高训练效率。

粗粒化

马尔科夫链不只是一个马尔科夫矩阵那么简单，还涉及到了不同的状态，也会有状态之间的转移概率，所以我们可以从不同的角度来看待粗粒化：

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

这三种出发点看似互有联系，但是具体关心的东西不同。第一种关注状态之间的连边拓扑结构，第二种关注不同状态的历史分布和相似性，第三种关注如何降维表达整个系统中的事件概率分布。

对马尔科夫链做粗粒化并没有一个统一的定义，目前该题目也没有一个完整的文献综述。在许多文献中，粗粒化coarse-graining，聚类partitioning/lumping/clustering和降维dimension reduction是重叠等价的。

多数文献中会出现如图一的定义：

图1：粗粒化定义

对状态空间[math]\displaystyle{ S }[/math]和转移函数[math]\displaystyle{ T(t):S \rightarrow S }[/math]，我们需要一个投影矩阵[math]\displaystyle{ R:S \rightarrow S' }[/math]，和一个新的转移函数[math]\displaystyle{ T'(t):S' \rightarrow S' }[/math]，而且[math]\displaystyle{ S' }[/math]和[math]\displaystyle{ T' }[/math]要符合马尔科夫链定义。

所以，我们能看出[math]\displaystyle{ R }[/math]是一个对状态空间[math]\displaystyle{ S }[/math]的降维函数（得到[math]\displaystyle{ S' }[/math]）。而且，我们还需要把[math]\displaystyle{ T }[/math]降维成[math]\displaystyle{ T' }[/math]。如果[math]\displaystyle{ T:s_{t+1} = s_{t} P }[/math]，[math]\displaystyle{ P }[/math]为离散马尔科夫矩阵的话，这个针对马尔科夫矩阵的降维可以表示为

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

(1)

其中[math]\displaystyle{ P \in \mathbb{R}^{n \times n} }[/math]，[math]\displaystyle{ P' \in \mathbb{R}^{k \times k} }[/math]，[math]\displaystyle{ W \in \mathbb{R}^{k \times n} }[/math]，[math]\displaystyle{ V \in \mathbb{R}^{n \times k} }[/math]。Buchholz的文章中把[math]\displaystyle{ W }[/math]称作distributor matrix，而[math]\displaystyle{ V }[/math]称作collector matrix。

马尔科夫链的粗粒化可以分成Hard partitioning和Soft partitioning。Hard partitioning是指把若干个微观状态分成一个宏观状态类，且一个微观状态不能同时属于多个宏观状态类，而Soft partitioning则会有可能出现一个微观状态被分在属于多个宏观类的情况。

在Hard partitioning的情况下，[math]\displaystyle{ V_{i, j} = 1\ if\ s_i \in s_j'\ and\ 0\ otherwise }[/math]，即一个表示[math]\displaystyle{ i }[/math]属于[math]\displaystyle{ j }[/math]列的零一分组矩阵，如果是Soft partitioning的话，[math]\displaystyle{ V }[/math]则不一定会是零一矩阵。

而

[math]\displaystyle{ \begin{aligned} W = D^{-1}(V)^T, D = diag(\alpha V), \end{aligned} }[/math]

(2)

其中[math]\displaystyle{ \alpha }[/math]是一个非负对角矩阵。

所以，我们能看到，对马尔科夫链做粗粒化的时候，不仅要对状态空间做降维([math]\displaystyle{ R }[/math])，还要对马尔科夫矩阵都要做降维([math]\displaystyle{ V }[/math])，当中也包含了对概率空间的降维([math]\displaystyle{ W }[/math])。

而且，有资料中还提到了粗粒化的交换律commutative。按照上面的定义，[math]\displaystyle{ s^{(t)} }[/math]有两个路径能达到[math]\displaystyle{ s'^{(t+1)} }[/math]：

1. 先得到下一步的[math]\displaystyle{ s^{(t+1)} = s^{(t)} P }[/math]，再做粗粒化得到[math]\displaystyle{ s'^{(t+1)} = R(s^{(t+1)}) }[/math]
2. 先做粗粒化得到[math]\displaystyle{ s'^{(t)} = R(s^{(t)}) }[/math]，再得到下一步的[math]\displaystyle{ s'^{(t+1)} = s'^{(t)} P' }[/math]

而交换律则要求从这两条路径得到的结果应该是一样的，即[math]\displaystyle{ R(s_{t} P) = R(s_{t}) P' }[/math]。

因此，根据目前有迹可循的对马尔科夫链的粗粒化的定义来说，我们要求一个良定义的粗粒化方法上述提到过的两个规则：

1. 粗粒化后的[math]\displaystyle{ S' }[/math]和[math]\displaystyle{ T' }[/math]要符合马尔科夫链定义。
2. 满足交换律

但是注意，这是我们从不同来源中总结出来的规则。目前我们还没找到一个严格的完备的定义。

而且，这只是最简单，最基本的对马尔科夫链的粗粒化，我们下面会介绍不同的指标，在上述的两个规则之上添加更多的规则会得出基于这些指标的粗粒化的方法。

在这个词条里，我们主要讨论Hard partitioning，主要参考的是Anru Zhang和Mengdi Wang的Spectral State Compression of Markov Processes，对Soft state aggregation感兴趣的朋友可以自行参阅S. P. Singh, T. Jaakkola, and M. I. Jordan, “Reinforcement learning with soft state aggregation."

Rank

首先讨论的是一个马尔科夫矩阵的rank秩。

大家理解的线代里的rank秩的定义是看矩阵中的线性无关的行向量的数量，但是这里对秩的理解是从一种类似于信道的概念。

秩的定义为我们能找到的一组概率密度函数 [math]\displaystyle{ f_1,\ ...\ ,\ f_r,\ g_1,\ ...\ ,\ g_r }[/math]，使得r在下列公式里最小：

[math]\displaystyle{ P(X_{t+1} | X_{t}) = \sum^r_{k=1} f_k(X_t) g_k(X_{t+1}) }[/math]

这里的秩的意思是，我们能多大程度上压缩信道，使得信息在宽度为秩的信道中无损传递。（笔者个人理解）

在n个离散状态的马尔科夫矩阵中，[math]\displaystyle{ f_1, ... , f_r, g_1, ... , g_r }[/math] 是维度为n的矩阵。

而我们能定义r × r的markov kernel [math]\displaystyle{ C = \{C_{ij} = \sum_{k=1}^n f_j(k)g_i(k)\} }[/math]，[math]\displaystyle{ f_1, ... , f_r }[/math] 为 left Markov features，[math]\displaystyle{ g_1, . . . , g_r }[/math] 为 right Markov features.

这个定义可以想象成可压缩的程度，也会是下面的hard partitioning的分组的数量。

Lumpability

Lumpability是一种描述状态分类的衡量，笔者暂时还没找到一个正式的中文翻译。

这个概念最早出现在Kemeny, Snell在1969年的Finite Markov Chains中。

初始定义

为了严谨的介绍这一定义，我们会详细翻译文章中的部分。不感兴趣的读者可以跳到更简洁直观的‘lumpable partition的充分必要条件'部分。

首先定义[math]\displaystyle{ s^{(t)} }[/math]表示系统在[math]\displaystyle{ t }[/math]时刻的微观状态，微观状态空间为[math]\displaystyle{ S=\{s_1, s_2, ... ,s_n\} }[/math]。

给定一个任意state partition [math]\displaystyle{ A=\{A_1, A_2, ... ,A_r\} }[/math]，也可以把其理解为宏观的状态空间。

[math]\displaystyle{ S }[/math]和[math]\displaystyle{ A }[/math]之间的Hard Partition映射关系为：[math]\displaystyle{ A_i \in S, A_i \neq \empty, A_i \cap A_j = \empty , \forall i, j, \cup_i A_i = S }[/math]

对于任意state partition [math]\displaystyle{ A }[/math]，我们能够定义一个lumped process，即把微观的动力学轨迹[math]\displaystyle{ s^{(t)} }[/math]投影到[math]\displaystyle{ A }[/math]的空间上。

这种轨迹的投影可以写作下列公式：

[math]\displaystyle{ \begin{aligned} &Pr_{\pi}[s^{(0)} \in A_i] \\ &Pr_{\pi}[s^{(1)} \in A_j | s^{(0)} \in A_i] \\ &Pr_{\pi}[s^{(t)} \in A_m | s^{(t-1)} \in A_k, ... , s^{(1)} \in A_j, s^{(0)} \in A_i] \end{aligned} }[/math]

(3)

其中，[math]\displaystyle{ \pi }[/math]为初始微观状态。

这些公式描述了：

1. 式子3.1：系统在[math]\displaystyle{ t=0 }[/math]的微观状态[math]\displaystyle{ s^{(0)} }[/math]属于某个宏观状态[math]\displaystyle{ A_i }[/math]的概率；
2. 式子3.2：已知在[math]\displaystyle{ t=0 }[/math]的微观状态[math]\displaystyle{ s^{(0)} }[/math]属于某个宏观状态[math]\displaystyle{ A_i }[/math]，那么在时间[math]\displaystyle{ t=1 }[/math]的微观状态[math]\displaystyle{ s^{(1)} }[/math]属于某个宏观状态[math]\displaystyle{ A_j }[/math]的概率；
3. 式子3.3：已知历史微观状态[math]\displaystyle{ \{s^{(0)}, ... , s^{(t-1)} \} }[/math]对应的宏观状态，那么在时间[math]\displaystyle{ t }[/math]的微观状态[math]\displaystyle{ s^{(t)} }[/math]属于某个宏观状态[math]\displaystyle{ A_m }[/math]的概率；

回想一下图(1)的交换律。我们现在有了微观动力学，也有了微观到宏观的粗粒化过程。所以式子(3)描述的是 微观状态->微观动力学->粗粒化->宏观状态 这一条路径。

需要注意的是，虽然任意一个state partition [math]\displaystyle{ A }[/math]都能走这样的路径，并形成lumped process。但不是所有的这些lumped process都是有意义且满足交换律的。

书中接下来就定义了什么叫做lumpable partition：

对于一个给定的state partition [math]\displaystyle{ A }[/math]，当

1. 式子(3)对任何微观初始状态(starting vector) [math]\displaystyle{ \pi }[/math] 都保持一致时，
2. 定义[math]\displaystyle{ A^{(t)} = A(s^{(t)}) }[/math]表示系统在[math]\displaystyle{ t }[/math]时刻的宏观状态，[math]\displaystyle{ \{A^{(0)}, ... , A^{(t-1)} \} }[/math]具有马尔科夫性，即[math]\displaystyle{ Pr_{\pi}[s^{(t)} \in A_m | s^{(t-1)} \in A_k, ... , s^{(1)} \in A_j, s^{(0)} \in A_i] = Pr_{\pi}[s^{(t)} \in A_m | s^{(t-1)} \in A_k] }[/math]时，

[math]\displaystyle{ A }[/math]是一个lumpable partition。

不满足上述两个条件的，都不被定义为lumpable partition。也就是说，不是所有的lumped process都是lumpable的，即使他们的命名方式相似。

lumpable partition的充分必要条件

对于lumpable partition，我们就能定义宏观动力学了。

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

设[math]\displaystyle{ p_{s_k \rightarrow s_m} = p(s^{(t)} = s_m | s^{(t-1)} = s_k) }[/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]\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^{(0)}, ... , A^{(t-1)} \} }[/math]，[math]\displaystyle{ p_{A_i \rightarrow A_j} = p(A^{(t)} = A_j | A^{(t-1)} = A_i) }[/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]中的所有状态的转移概率的和。

对于lumpable partition [math]\displaystyle{ A }[/math]定义宏观动力学后，我们就能走另外一条路径了：微观状态->粗粒化->宏观动力学->宏观状态。

请注意：lumpability有一个非常重要的前提，即我们要给定一个对马尔科夫链的state的partition，然后根据马尔科夫状态转移矩阵来决定这个partition对这个马尔科夫链是否lumpable。

所以，lumpability并不是一个马尔科夫链的特质，而是对马尔科夫链的某个state partition的特质。我们不能直接说某个马尔科夫链是lumpable或non-lumpable，但是我们可以说，对于某个马尔科夫链的某几个状态partition是lumpable或non-lumpable。

Lumpability正例和反例

对于一个马尔科夫链，有些分组partition是lumpable，有些分组是non-lumpable。这里用一个简单的例子提供正例和反例。

我们用到的是这样一个转移矩阵：

[math]\displaystyle{ \begin{aligned} P = P_{s_i \rightarrow s_j} = \left [ \begin{array}{ccc:c} 0.3 & 0.3 & 0.3 & 0.1 \\ 0.3 & 0.3 & 0.3 & 0.1 \\ 0.3 & 0.3 & 0.3 & 0.1 \\ \hdashline0.1 & 0.1 & 0.1 & 0.7 \end{array} \right ] \end{aligned} }[/math]

正例

我们可以看到，123行都是一样的，所以把他们[math]\displaystyle{ \{s_1, s_2, s_3\} }[/math]分在一组[math]\displaystyle{ A_1 }[/math]，并把状态[math]\displaystyle{ s_4 }[/math]分在另外一组[math]\displaystyle{ A_2 }[/math]，看起来像是一个合理的正例。

接下来我们来计算一下。

设[math]\displaystyle{ p_{s_1 \rightarrow s_2} = p(s^{(t)} = s_2 | s^{(t-1)} = s_1) }[/math]，[math]\displaystyle{ p_{s_1 \rightarrow A_2} = p(s^{(t)} \in A_2 | s^{(t-1)} = s_1) }[/math]

第1组[math]\displaystyle{ A_1 }[/math]的状态[math]\displaystyle{ s_1 }[/math]到第1组[math]\displaystyle{ A_1 }[/math]，[math]\displaystyle{ p_{s_1 \rightarrow A_1} = p_{A_1 \rightarrow A_1} = \sum_{s_m \in A_1} p_{s_1 \rightarrow s_m} = p_{s_1 \rightarrow s_1} + p_{s_1 \rightarrow s_2} + p_{s_1 \rightarrow s_3} = 0.3 \times 3 = 0.9 }[/math]

第1组[math]\displaystyle{ A_1 }[/math]的状态[math]\displaystyle{ s_1 }[/math]到第2组[math]\displaystyle{ A_2 }[/math]，[math]\displaystyle{ p_{s_1 \rightarrow A_2} = p_{A_1 \rightarrow A_2} = \sum_{s_m \in A_2} p_{s_1 \rightarrow s_m} = p_{s_1 \rightarrow s_4} = 0.1 }[/math]

第2组[math]\displaystyle{ A_2 }[/math]的状态[math]\displaystyle{ s_4 }[/math]到第1组[math]\displaystyle{ A_1 }[/math]，[math]\displaystyle{ p_{s_4 \rightarrow A_1} = p_{A_2 \rightarrow A_1} = \sum_{s_m \in A_1} p_{s_4 \rightarrow s_m} = p_{s_4 \rightarrow s_1} + p_{s_4 \rightarrow s_2} + p_{s_4 \rightarrow s_3} = 0.1 \times 3 = 0.3 }[/math]

第2组[math]\displaystyle{ A_2 }[/math]的状态[math]\displaystyle{ s_4 }[/math]到第2组[math]\displaystyle{ A_2 }[/math]，[math]\displaystyle{ p_{s_4 \rightarrow A_2} = p_{A_2 \rightarrow A_2} = \sum_{s_m \in A_2} p_{s_4 \rightarrow s_m} = p_{s_4 \rightarrow s_4} = 0.7 }[/math]

[math]\displaystyle{ \begin{aligned} P = P_{A_k \rightarrow A_l} = \left [ \begin{array}{c:c} 0.9 & 0.1 \\ \hdashline0.3 & 0.7 \end{array} \right ] \end{aligned} }[/math]

反例

现在我们来看一个反例：把状态12分为一组，状态34分为另一组，

同样的，我们来计算一下。这次我们需要对3和4分开计算了。

第1组[math]\displaystyle{ A_1 }[/math]的状态[math]\displaystyle{ s_1 }[/math]到第1组[math]\displaystyle{ A_1 }[/math]，[math]\displaystyle{ p_{s_1 \rightarrow A_1} = p_{A_1 \rightarrow A_1} = \sum_{s_m \in A_1} p_{s_1 \rightarrow s_m} = p_{s_1 \rightarrow s_1} + p_{s_1 \rightarrow s_2} = 0.3 \times 2 = 0.6 }[/math]

第1组[math]\displaystyle{ A_1 }[/math]的状态[math]\displaystyle{ s_1 }[/math]到第2组[math]\displaystyle{ A_2 }[/math]，[math]\displaystyle{ p_{s_1 \rightarrow A_2} = p_{A_1 \rightarrow A_2} = \sum_{s_m \in A_2} p_{s_1 \rightarrow s_m} = p_{s_1 \rightarrow s_3} + p_{s_1 \rightarrow s_4} = 0.4 }[/math]

第2组[math]\displaystyle{ A_2 }[/math]的状态[math]\displaystyle{ s_3 }[/math]到第1组[math]\displaystyle{ A_1 }[/math]，[math]\displaystyle{ p_{s_3 \rightarrow A_1} = \sum_{s_m \in A_1} p_{s_3 \rightarrow s_m} = p_{s_3 \rightarrow s_1} + p_{s_3 \rightarrow s_2} = 0.3 \times 2 = 0.6 }[/math]

第2组[math]\displaystyle{ A_2 }[/math]的状态[math]\displaystyle{ s_4 }[/math]到第1组[math]\displaystyle{ A_1 }[/math]，[math]\displaystyle{ p_{s_4 \rightarrow A_1} = \sum_{s_m \in A_1} p_{s_4 \rightarrow s_m} = p_{s_4 \rightarrow s_1} + p_{s_4 \rightarrow s_2} = 0.1 \times 2 = 0.2 }[/math]

这里我们就能看到，[math]\displaystyle{ p_{s_3 \rightarrow A_1} \neq p_{s_4 \rightarrow A_1} }[/math]，当同一组的两个状态[math]\displaystyle{ s_3 }[/math]和[math]\displaystyle{ s_4 }[/math]对其他组的转移概率不一样的话，如果我们强行按照这样来分组（也没办法强行，因为我们不知道[math]\displaystyle{ p_{A_2 \rightarrow A_1} = p_{s_3 \rightarrow A_1} }[/math]还是[math]\displaystyle{ p_{A_2 \rightarrow A_1} = p_{s_4 \rightarrow A_1} }[/math]），假设[math]\displaystyle{ p_{A_2 \rightarrow A_1} = 0.6 }[/math]，我们会发现这样的粗粒化结果违背了Lumpability一开始的定义（公式(3)），即粗粒化后的转移概率对所有的初始微观状态[math]\displaystyle{ \pi }[/math]都适用。因为当[math]\displaystyle{ \pi = s_4 }[/math]的时候，[math]\displaystyle{ p_{A_2 \rightarrow A_1} = 0.6 }[/math]这个转移概率就是错的，即使[math]\displaystyle{ \pi \neq s_4 }[/math]，任何初始微观状态[math]\displaystyle{ \pi }[/math]都会在某时刻[math]\displaystyle{ n }[/math]达到[math]\displaystyle{ s_4 }[/math]并导致转移概率出错。

一个更复杂的正例

三个相同状态的例子可能稍微有点简单，接下来我们来看一个稍微复杂一点的例子：

[math]\displaystyle{ \begin{aligned} P = P_{s_i \rightarrow s_j} = \left [ \begin{array}{cc:cc} 0.2 & 0.4 & 0.3 & 0.1 \\ 0.4 & 0.2 & 0.1 & 0.3 \\ \hdashline0.2 & 0.1 & 0.2 & 0.5 \\ 0.1 & 0.2 & 0 & 0.7 \end{array} \right ] \end{aligned} }[/math]

对于这个马尔科夫矩阵，[math]\displaystyle{ \{\{1,2\}, \{3,4\}\} }[/math]是一个lumpable的分组。分组后的转移矩阵为： [math]\displaystyle{ \begin{aligned} P = P_{A_k \rightarrow A_l} = \left [ \begin{array}{c:c} 0.6 & 0.4 \\ \hdashline0.3 & 0.7 \end{array} \right ] \end{aligned} }[/math]

Lumpability和粗粒化

我们在粗粒化的部分提到了，马尔科夫链的粗粒化不仅要对状态空间做，也要对转移矩阵和概率空间做。

这三个部分可以同时做，也可以看作为：先对状态空间做，再对转移矩阵和概率空间做；即先对状态做分组，然后再获取对应的粗粒化后的转移矩阵。

我们也提到过对状态空间做粗粒化有Hard Partitioning和Soft Partitioning两种。Soft Partitioning可以看作把微观状态打散重构成了一些宏观状态，而Hard Partitioning则是更严格的，把若干个微观状态分成一个组。

而Lumpability就是一个指标，用来评价‘对于某一种Hard Partitioning的微观状态分组分案，是否对微观状态转移矩阵lumpable’。

不管状态空间按照哪一个Hard Partitioning方案做分类，它都有对应后续的对转移矩阵和概率空间的粗粒化方案（公式(1)和(2)），其中有些分组满足交换律，有些则不满足。

其中的某些分组方案lumpable，也有某些分组方案non-lumpable。

所以，整个链条应该是这样的：

所以，一个满足lumpability的良定义的粗粒化方案应该在原来的两个条件下，再多加两个条件：

1. 粗粒化后的[math]\displaystyle{ S' }[/math]和[math]\displaystyle{ T' }[/math]要符合马尔科夫链定义；
2. 满足交换律；
3. Hard Partitioning，即存在分组矩阵[math]\displaystyle{ V_{i, j} = 1\ if\ s_i \in s_j'\ and\ 0\ otherwise }[/math]；
4. [math]\displaystyle{ V }[/math]是一个lumpable分组。

Lumpable Partition粗粒化的一致性Consistency

这一小节里我们讨论一下，对一个马尔科夫链的lumpable partition 与 HOM（Higher-order Macronodes，详情请参考复杂网络中的因果涌现），和Consistency之间的关联。

Higher-order Macronodes，是由Jian Xu提出，根据给定的微观马尔科夫链和节点的聚类方案，得出一个和微观转移矩阵保持一致的宏观转移矩阵的方法。

我们这里使用一个简单的lumpable partition例子，来显示lumpability的粗粒化跟HON是对应的。

图2：例子示意图

[math]\displaystyle{ \begin{aligned} P = P_{s_i \rightarrow s_j} = \left [ \begin{array}{c:c:cc:c:c} 0 & 0 & 0.3 & 0.1 & 0.6 & 0\\ \hdashline 0 & 0 & 0.2 & 0.3 & 0 & 0.5\\ \hdashline 0 & 0 & 0 & 0.7 & 0.2 & 0.1\\ 0 & 0 & 0.7 & 0 & 0.2 & 0.1\\ \hdashline 1.0 & 0 & 0 & 0 & 0 & 0\\ \hdashline 0 & 1.0 & 0 & 0 & 0 & 0\\ \end{array} \right ] \end{aligned} }[/math]

对于这个转移概率矩阵，[math]\displaystyle{ \{1, 2, \{34\}, 5, 6\} }[/math]是一个lumpable partition。

首先，我们先按照lumpability的粗粒化定义，计算对应的宏观转移概率： [math]\displaystyle{ \begin{aligned} \left [ \begin{array}{c:c:c:c:c} 0 & 0 & 0.4 & 0.6 & 0\\ \hdashline 0 & 0 & 0.5 & 0 & 0.5\\ \hdashline 0 & 0 & 0.7 & 0.2 & 0.1\\ \hdashline 1.0 & 0 & 0 & 0 & 0\\ \hdashline 0 & 1.0 & 0 & 0 & 0\\ \end{array} \right ] \end{aligned} }[/math]

然后，我们来简单介绍一下的构建宏观HOM转移概率矩阵的计算方法。

在图2中，我们看到3和4为待合并节点[math]\displaystyle{ A }[/math]。

1. 对于其他节点到待合并节点[math]\displaystyle{ A }[/math]的连边，即待合并节点的入流（图2中绿色的连边），我们直接加和即可，即[math]\displaystyle{ p_{s_1 \rightarrow A} = 0.3+0.1=0.4 }[/math]，[math]\displaystyle{ p_{ s_2 \rightarrow A} = 0.2+0.3=0.5 }[/math]

2. 对于待合并节点[math]\displaystyle{ A }[/math]到其他节点的连边，即待合并节点的出流（图2中棕色的连边），我们需要考虑三种情况：

1. 当待合并节点[math]\displaystyle{ A }[/math]间存在连边时；
2. 当待合并节点[math]\displaystyle{ A }[/math]指向同一个输出节点时；
3. 当待合并节点[math]\displaystyle{ A }[/math]指向不同输出节点时；

待合并节点[math]\displaystyle{ A }[/math]的出流的计算方式都不一样，如图3所示。

但是，我们可以注意到，这三种计算方式都是待合并节点[math]\displaystyle{ A }[/math]中的各节点的出流的[math]\displaystyle{ W_i^{out} }[/math]加权，而我们知道，lumpability的特性决定了，群组里的各节点到其他群组的出流[math]\displaystyle{ W_i^{out} }[/math]都是一样的，对相同的[math]\displaystyle{ W_i^{out} }[/math]做任意加权平均的结果都等于[math]\displaystyle{ W_i^{out} }[/math]。也就是说，任意的加权方法对于lumpable partition来说都能得出相同的结果。

由此，我们总结，lumpable partition的聚合方式和HOM的计算方式，无论入流还是出流都完全一致。

我们得出这个结论的目标是，借助HOM方法的consistency的性质，说明lumpable partition的聚合方式，也有consistency的性质。

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

定义微观网络[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），若结果为零则动力学一致。公式为：

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

也就是说，Inconsistency 等于 宏观和微观网络节点在不同时刻的概率分布的KL散度之和。

实验发现，对不同节点规模以及参数下的微观网络，使用HOM粗粒化后的宏观网络的不一致性随着迭代步数的增加都会收敛到0。这说明，HOM构建的宏观网络和微观网络是保持一致的。

回到上面lumpable partition 等于HOM的推断，我们就能得出，根据lumpability而做的粗粒化，构建的宏观网络和微观网络也是保持一致的。

基于Lumpability的粗粒化方法（未给定lumpable partition的情况）

图4：Zhang 文章中的示意图。图中左面四个矩阵都是lumpable马尔科夫矩阵，而右面的P_2是一个噪声矩阵，(P_1)^T P_2 = 0

由上面的lumpability公式(4)和例子中我们能获得一个直观上的说法：当马尔科夫矩阵存在block结构，或者状态明显可被分成几类的时候，根据这样的partition，该矩阵就会lumpable，如图4中的[math]\displaystyle{ \bar{P} }[/math]所示，把相同的状态（行向量）分成一类的partition显然lumpable。

但是，有时候有些lumpable的矩阵的状态排序被打乱了（如图一中的[math]\displaystyle{ P_1 }[/math]），或者矩阵包含了如[math]\displaystyle{ P_2 }[/math]的噪声（如图4中的[math]\displaystyle{ P }[/math]，[math]\displaystyle{ P = P_1 + P_2 }[/math]，[math]\displaystyle{ P_1^TP_2 = 0 }[/math]）。

我们在实际问题中很多时候要面对的是像[math]\displaystyle{ P }[/math]这样的矩阵，我们既无法确定它是否lumpable，也无法决定它的partition，我们甚至不知道它的马尔科夫秩。

在这种情况下，Anru Zhang的文章中提供了一种寻找最优partition的方法。通过此方法我们能找到最优的partition，代入这个partition我们就能通过上面的充分必要条件来决定一个马尔科夫矩阵是否lumpable。具体步骤如下：

1. 先获取一个n*n维的马尔科夫矩阵P，或者是马尔科夫链的频率采样；
2. 获取r<n 的马尔科夫秩，或指定一个r；
3. 对P进行SVD分解，[math]\displaystyle{ P = U \Sigma V^T }[/math]，其中U为左奇异向量，V为右奇异向量。
4. 通过下列公式得到最优partition [math]\displaystyle{ \Omega_1, \Omega_2, ... , \Omega_r }[/math]

[math]\displaystyle{ \Omega_1, \Omega_2, ... , \Omega_r = argmin_{\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]类[math]\displaystyle{ \Omega_s }[/math]的特征向量。

这个算法的意思是在最优的partition中，[math]\displaystyle{ \Omega_s }[/math]中的状态[math]\displaystyle{ i }[/math]的右奇异向量会和[math]\displaystyle{ v_s }[/math]距离最小，也就是说，让每个微观态的右特征向量和它对应的宏观态的右特征向量尽可能接近。

这个看起来非常眼熟，其实它就是在做kMeans聚类算法。让[math]\displaystyle{ v_s }[/math]和组里的[math]\displaystyle{ V[i:r] }[/math]距离最小的方法，就是让[math]\displaystyle{ v_s }[/math]成为[math]\displaystyle{ V[i:r] }[/math]这若干个点的中心。

回顾一下kMeans算法把n个点聚成r类的目标函数，其中第[math]\displaystyle{ i }[/math]个点被归到第[math]\displaystyle{ c^i }[/math]类，而[math]\displaystyle{ \mu_j }[/math]为第[math]\displaystyle{ j }[/math]类点的中心点：

[math]\displaystyle{ J(c^1, ... , c^n, \mu_1, ... , \mu_r) = \sum_{i=1}^n || x_i - \mu_{c^i} ||^2 }[/math]

我们看到这两个目标函数是非常相似的。所以，我们就可以通过上述算法或者kMeans来获取最优partition，然后根据这个partition来判断马尔科夫矩阵的lumpability。

因果涌现

有效信息

有效信息相关详情请参照因果涌现词条。

Dynamical Reversibility 动力学可逆性

动力学可逆性相关详情请参照因果涌现词条。

因果态

因果态相关详情请参照计算力学词条。

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

动力学可逆性 v.s. Lumpability

在lumpability的图1中我们提到了lumpable的马尔科夫矩阵可以被重新排列成几个block，这种lumpable的矩阵的动力学可逆性也会很高，在这种情况下动力学可逆性和 Lumpability是一致的。

在没有明显block结构的情况下，我们可以使用上述的SVD+kMeans方式找到最优的lumpability partition并决定它的lumpability。我们在一个动力学可逆性比较高（有明显的因果涌现）的国际贸易网的例子上尝试了一下，却发现该网络上的最优lumpability partition宏观结果的动力学可逆性并不高。这现象说明：

1. Lumpability最优的partition不等于可逆性最优的partition。

1. Lumpability的分组数量是基于马尔科夫秩的，而可逆性的分组数量是对奇异值谱的截断位置。

（未完待续）