基于马尔科夫链的流网络分析

一个平衡的流网络可以自然地转变成马尔科夫链，而马尔科夫链是数学家们研究得非常深入而透彻的数学对象，因此这为我们更好地分析流网络提供了便利条件。 下面给出基本计算简表。定义：流矩阵F，对应的马尔科夫链矩阵M [math]\displaystyle{ m_{ij}=\frac{f_{ij}}{\sum_{j=1}^{N+1}f_{ij}} }[/math] 基础矩阵:[math]\displaystyle{ U=I+M+M^2+\cdot\cdot\cdot=(I-M)^{-1} }[/math]

[math]\displaystyle{ m_{ij}=\frac{f_{ij}}{\sum_{j=1}^{N+1}f_{ij}} }[/math]，

基础矩阵[math]\displaystyle{ U=I+M+M^2+\cdot\cdot\cdot=(I-M)^{-1} }[/math] i到j的总流：[math]\displaystyle{ t_{i,j}=T_0\frac{u_{0,i}u_{i,j}}{u_{i,i}} }[/math]，i到j的首达流：[math]\displaystyle{ \phi_{i,j}=T_0 \frac{u_{0,i}u_{i,j}}{u_{i,i}u_{j,j}} }[/math]

其中T₀为从源输入到整个网络的流量，也就是IS（参看流动网络）。根据这两个公式，可以计算任意网络任意节点对的总流和首达流。

i到j的平均时间：[math]\displaystyle{ k_{i,j}=\frac{(MU^2)_{ij}}{(MU)_{ij}} }[/math]，i到j的首达平均时间：[math]\displaystyle{ l_{i,j}=\frac{u_{jj}(M_{-j}U_{-j}^2)_{ij}}{u_{ij}}=\frac{(MU^2)_{ij}}{u_{ij}}-\frac{(MU^2)_{jj}}{u_{jj}} }[/math]

其中[math]\displaystyle{ M_{-j} }[/math]表示将M的第i行全置为0，[math]\displaystyle{ U_{-j}=I+M_{-j}+M_{-j}^{2}+\cdot\cdot\cdot=(I-M)^{-1} }[/math]。

从流网络到马尔科夫链的转换

图1：一个示例流网络

对于平衡流动来说，我们总可以将整个流系统看作是一个多粒子沿着网络流动的系统。例如，对于图1所示网络，当某个处于1号节点的粒子往外流动的时候，它可能会跳到2号节点，也可能会跳到3号节点。可以想到，由于1到2号节点的流量比到3号的节点大，所以这个粒子更有可能跳到2号节点。我们不妨假设粒子跳转到某个节点的概率是与这条边上的流量成正比的。也就是说，粒子会以5/8=50/80的概率从1跳到2，而以3/8=30/80的概率从1到3。同样的道理，当粒子跳入到任何一个节点后，它都会继续沿着边以概率跳转。

我们可以写出这个网络所对应的流矩阵（参见流动网络）如下：

其中，矩阵外围的数字为对节点的编码。0对应源，6对应汇。阴影区域对应网络的实线部分，为我们关注的节点主体。 整个平衡态的流网络可以看作是一个马尔可夫链，其中每个节点相当于是粒子所处的可能状态，节点之间的跳转可以看作是概率转移。具体的，对于任何的处于平衡的流网络F，我们定义矩阵M为概率转移矩阵，其中M中的第i行第j列的元素就是粒子处于i状态下转移到j状态的概率，记为：

[math]\displaystyle{ m_{ij}=Pr\{to~~j|after~~visit~~i\}=\frac{f_{ij}}{\sum_{j=1}^{N+1}f_{ij}} }[/math]

当粒子进入汇节点N+1，它就不会再往外跳了，也就是说汇所在的行全为0，因此，状态N+1就相当于马尔可夫链的吸收态。注意，根据矩阵M的定义，它自然满足概率归一化条件（不包括汇节点）：

[math]\displaystyle{ \sum_{j=1}^{N+1}m_{ij}=\frac{\sum_{j=1}^{N+1}f_{ij}}{\sum_{j=1}^{N+1}f_{ij}}=1,~~~~\forall i\lt N+1 }[/math]

针对图1的示例网络，转化过来的马尔科夫矩阵为：

其中阴影部分为我们主要关心的马尔科夫转移概率。

逆向马尔可夫链

另外一个比较有趣的事实是，根据流矩阵F，我们还可以定义另外一个马尔科夫链M'，称为逆向马尔科夫链。它刻画的是一个逆向的因果过程，也就是说如果我在图1中2号节点观察到了一个新到的粒子，那么它有可能从1号节点过来也可能从3号节点过来。显然1到2的流量要比3到2的流量大得多，因此粒子从1跳过来的可能性会更大。那么，既然整个网络是平稳的，我们不难想到，在2号节点观测到一个粒子的条件下，它是从1号节点跳过来的概率是5/6=50/60，它从3号节点跳过来的概率是1/6=10/60。

更一般地，我们定义逆向马尔科夫链M'，其中任意的元素m'_ij定义为：观察到粒子来到了节点i的条件下，它可能是从j节点流过来的概率。它按照下式计算：

[math]\displaystyle{ m'_{ij}=Pr\{from~~j|observed~~particle~~on~~i\}=\frac{f_{ji}}{\sum_{j=0}^{N}f_{ji}} }[/math]

同样地道理，它也是归一化的，即：

[math]\displaystyle{ \sum_{j=0}^{N}m'_{ji}=1,~~~~\forall i\gt 0 }[/math]

这个逆向马尔可夫链又可以看作是逆向流矩阵F^T(即原始流矩阵F的转置)所对应的正向马尔科夫矩阵，这是因为：

[math]\displaystyle{ m'_{ij}=\frac{f_{ji}}{\sum_{j=0}^{N}f_{ji}}=\frac{F^T_{\{ij\}}}{\sum_{j=0}^{N}F^T_{\{ij\}}}=m_{ij}(F^T) }[/math]

其中逆向流矩阵就相当于把原始流网络所有的有向边都掉转方向，并保持所有的流量不变。由于网络是平衡的，所以逆向的流网络必然也是平衡的。

这似乎给我们描绘了一个有趣的图景：在一个平衡的流网络中，所有的粒子都顺着流动从源到汇。而与此同时，如果我们从汇观察到粒子，并不断地询问这个粒子是从哪个节点跳转过来的，我们就得到了一个反向的信息流动，就仿佛是有假想的粒子从汇流到了源一样。

路径分析

下面，我们沿着粒子在流网络中流动，分析粒子所能经过的所有可能路径。让我们以图1所示网络为例，考察从1节点到4节点的路径。我们可以将所有从1到4的路径分成2类：从1到4的首达路径，以及包含了从4到4的循环路径。而对于首达路径，又可以分成2类：直接路径、首达流循环路径。

首达路径

首先，我们可以列出从1到4的所有直接路径，所谓的直接路径是指粒子从1流到4，路径上的所有节点只被访问一次（没有循环）的所有可能路径，对于这个简单的网络来说，我们可以直接列出这些路径是：

[math]\displaystyle{ 1\rightarrow 2\rightarrow 4 }[/math]

[math]\displaystyle{ 1\rightarrow 3\rightarrow 2\rightarrow 4 }[/math]

其次，从1到4的首达流循环路径是指粒子从1流到4，中间允许任意次地循环重复经过某个节点，但是一旦到达4之后就不再循环重复的所有可能路径。也就是说，4节点只能唯一地出现在路径的尾端。对于图1所示的网络，我们知道1-->2-->5-->3-->2-->4是一条首达循环路径，而1-->2(-->5-->3-->2)(-->5-->3-->2)-->4也是一条首达循环路径，其中括号的部分被重复了两次。当然，括号部分中重复3次、4次…… 都是不同的首达循环路径。所以，我们可以把这一系列从5-->3-->2之间循环任意次的路径简记成

[math]\displaystyle{ 1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m\rightarrow 2\rightarrow 4 }[/math]

其中，指数的意思为括号部分的路径字符可以重复m次，m为任意的整数（可以取无穷大），因此有无穷条从1到4的首达循环路径。

值得指出的是，在这里，我们相当于定义了两个路径的乘法，如果路径P1=i1-->...-->i，路径P2=j1,....,j，则P1乘以P2定义为：

[math]\displaystyle{ P_1 \cdot P_2=i_1\rightarrow i_2 \rightarrow \cdot\cdot\cdot i \rightarrow j_1\rightarrow j_2\rightarrow \cdot\cdot\cdot j }[/math]

可以验证上述的路径平方就是我们定义的路径的乘法。

同样的道理，还有一系列不同的从1到4的首达循环路径：

综上所述，我们实际上可以把所有的从1到4的首达路径（包含了不循环的和循环的）统一表示成：

[math]\displaystyle{ 1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m\rightarrow 4,~~~~ 1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2\rightarrow 4 }[/math]

其中m为大于等于0的任意整数。当m=0的时候，就是直接路径。

循环路径

第二大类路径是那些到了4号节点又沿着网络的流动重复到达4号节点的流动路径，它们就是循环路径，例如：1-->2-->4-->5-->3-->2-->4，又如：1-->3-->2-->5-->3-->4-->5-->3-->2-->4。利用我们上面引入的记号，我们可以把这两个路径进行推广统一写为：

[math]\displaystyle{ [1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m\rightarrow 4]\cdot[(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^t }[/math]

[math]\displaystyle{ [1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2\rightarrow 4]\cdot [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^t }[/math]

其中m,m',t都是大于等于0的任意整数。我们看到，上式两种路径就是在所有可能的首达路径后面“后缀”了一个从4到4的首达路径的任意次方。如果我们再定义路径的加法为将两个相加的路径并列在一起（相当于集合的并集），则我们可以简洁统一地把从1到4的所有可能路径表达为：

[math]\displaystyle{ \sum_{m,m',t}{[1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m \rightarrow 4]\cdot [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^t+ [1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2\rightarrow 4]\cdot [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^{t}} }[/math]

可以验证路径的加法和乘法同样满足结合律，所以上式又可以写作：

[math]\displaystyle{ \sum_{m,m',t}{\{[1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m ] + [1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2]\}\cdot [(\rightarrow 4\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^{t}} }[/math]

(eq1)

对于这种简单的网络来说，要表达出所有从1到4的路径都已经如此复杂，当我们考虑更复杂的网络的时候，机械地列出所有路径就几乎是不可能的了。所以，我们必须找出一套更合理的方法。

路径矩阵

我们可以用矩阵来表示路径，这样，通过矩阵的乘积可以表示路径的合成。首先，根据原始网络的邻接矩阵，我们定义一个矩阵P：

[math]\displaystyle{ P_{ij}=\left\{\begin{array}{ll} i\rightarrow j & \mbox {if i,j connected}, \\ 0 &\mbox {else} \end{array}\right. }[/math]

也就是，如果i和j相邻，则在(i,j)的位置上放置i-->j，表示从i到j的路径。该矩阵就是路径矩阵。

例如，图1所示网络的路径矩阵可以写为：

[math]\displaystyle{ P_1=\begin{pmatrix} 0 &0\rightarrow 1&0 &0 &0 &0 &0 & \\ 0 &0 &1\rightarrow 2&1\rightarrow 3&0 &0 &0 & \\ 0 &0 &0 &0 &2\rightarrow 4&2\rightarrow 5&2\rightarrow 6& \\ 0 &0 &3\rightarrow 2&0 &0 &0 &3\rightarrow 6& \\ 0 &0 &0 &0 &0 &4\rightarrow 5&4\rightarrow 6& \\ 0 &0 &0 &5\rightarrow 3&0 &0 &5\rightarrow 6& \\ 0 &0 &0 &0 &0 &0 &0 & \\ \end{pmatrix} }[/math]

下面我们来考虑路径矩阵的乘积。设路径矩阵P和Q乘积（它们都是n*n的方阵）得到R，则

[math]\displaystyle{ R_{ij}=\sum_{k}P_{ik}\cdot Q_{kj} }[/math]

也就是说，路径矩阵的乘积与普通矩阵一样。注意等式右侧的乘法对应路径的乘法，也就是如果[math]\displaystyle{ P_{ik}=i\rightarrow k, Q_{kj}=k\rightarrow j }[/math]，则：[math]\displaystyle{ P_{ik}\cdot Q_{kj}=i\rightarrow k\rightarrow j }[/math]，即把两个路径拼接在一起。

这样P₁自乘一次就得到：

[math]\displaystyle{ \begin{array}{ll} P_1^2=\begin{pmatrix} 0 &0\rightarrow 1&0 &0 &0 &0 &0 & \\ 0 &0 &1\rightarrow 2&1\rightarrow 3&0 &0 &0 & \\ 0 &0 &0 &0 &2\rightarrow 4&2\rightarrow 5&2\rightarrow 6& \\ 0 &0 &3\rightarrow 2&0 &0 &0 &3\rightarrow 6& \\ 0 &0 &0 &0 &0 &4\rightarrow 5&4\rightarrow 6& \\ 0 &0 &0 &5\rightarrow 3&0 &0 &5\rightarrow 6& \\ 0 &0 &0 &0 &0 &0 &0 & \\ \end{pmatrix}\cdot \begin{pmatrix} 0 &0\rightarrow 1&0 &0 &0 &0 &0 & \\ 0 &0 &1\rightarrow 2&1\rightarrow 3&0 &0 &0 & \\ 0 &0 &0 &0 &2\rightarrow 4&2\rightarrow 5&2\rightarrow 6& \\ 0 &0 &3\rightarrow 2&0 &0 &0 &3\rightarrow 6& \\ 0 &0 &0 &0 &0 &4\rightarrow 5&4\rightarrow 6& \\ 0 &0 &0 &5\rightarrow 3&0 &0 &5\rightarrow 6& \\ 0 &0 &0 &0 &0 &0 &0 & \\ \end{pmatrix}&\\ = \begin{pmatrix} 0 &0 &0\rightarrow 1\rightarrow 2&0\rightarrow 1\rightarrow 3 &0 &0 &0 & \\ 0 &0 &0 &1\rightarrow 3\rightarrow 2&1\rightarrow 2\rightarrow 4&1\rightarrow 2\rightarrow 5&(1\rightarrow 2\rightarrow 6)+(1\rightarrow 3\rightarrow 6)& \\ 0 &0 &0 &2\rightarrow 5\rightarrow 3&0&2\rightarrow 4\rightarrow 5&(2\rightarrow 4\rightarrow 6)+(2\rightarrow 5\rightarrow 6)& \\ 0 &0 &0 &0 &3\rightarrow 2\rightarrow 4&3\rightarrow 2\rightarrow 5&3\rightarrow 2\rightarrow 6& \\ 0 &0 &0 &4\rightarrow 5\rightarrow 3&0 &0 &4\rightarrow 5\rightarrow 6& \\ 0 &0 &5\rightarrow 3\rightarrow 2&0 &0 &0 &5\rightarrow 3\rightarrow 6& \\ 0 &0 &0 &0 &0 &0 \\ \end{pmatrix} &\end{array} }[/math]

其中i,j元素表示从i到j长度为2的所有路径。同样的道理，我们可以计算P₁^m，它的任意元素就表示从i到j长度为m的所有可能路径。从这点来看，矩阵的乘积与路径的合成有着天然的对应。

全路径

接下来我们就要利用路径矩阵这个工具计算从任意节点i到任意节点j的所有可能路径。我们知道P^m的第i,j个元素包括了从i到j的长度为m的可能路径，那么要得到从i到j的所有可能路径（无论路径有多长），我们只要把所有的P¹、P²……P^m……加起来就可以了。所以：

[math]\displaystyle{ Paths\{i\rightarrow j\}=(\sum_{m=1}^{\infty}P^m)_{ij} }[/math]

例如，对于图1的网络，1到4的所有可能路径就是：

[math]\displaystyle{ Paths\{1\rightarrow 4\}=(\sum_{m=1}^{\infty}P_1^m)_{1,4} }[/math]

如果我们把P₁代入到上式，可以验证它跟我们在eq1中得到的全路径的表达式是一样的。

除此之外，我们也可以写出所有从i到j的首达路径（j在整个路径中只能出现在结尾）为：

[math]\displaystyle{ First\{i\rightarrow j\}=\frac{(\sum_{m=1}^{\infty}P^m)_{i,j}}{(\sum_{m=1}^{\infty}P^m)_{j,j}} }[/math]

这里，两个路径相除定义为乘法的逆运算。即，如果[math]\displaystyle{ R\cdot Q=P }[/math]，则定义[math]\displaystyle{ \frac{P}{Q}=R }[/math]。例如，若[math]\displaystyle{ P=1\rightarrow 2\rightarrow 3\rightarrow 4, Q=3\rightarrow 4 }[/math], 则[math]\displaystyle{ \frac{P}{Q}=1\rightarrow 2 }[/math]

我们可以写出图1中所有从1到4的首达路径为：

[math]\displaystyle{ First\{1\rightarrow 4\}=\frac{(\sum_{m=1}^{\infty}P_1^m)_{1,4}}{(\sum_{m=1}^{\infty}P_1^m)_{4,4}} }[/math]

根据eq1，我们知道所有的首达路径拼接上（乘上）任意一条从4到4的路径就构成了一条非首达路径，那么我们将所有的路径除以一条从4到4的路径就可以得到首达路径了。也就是在eq1中将[math]\displaystyle{ [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^{t} }[/math]这个因子除去。

流量分析

总流量

我们知道，每一条从i到j的路径就代表了粒子从i到j的一种可能流动。当有大量的粒子在体系中流动的时候，每条路径上都会分配一定的流量。下面，我们就要求出从i到j经过所有可能路径的粒子总流量[math]\displaystyle{ t_{i,j} }[/math]。这个流量可以通过这样的假想实验来定义：假设原始的粒子没有颜色，而粒子一旦经过节点i就被染成了红色，并且一直保持不褪色。那么，我从j节点处统计：每一时刻流入j节点的红色粒子数有多少？这个数量就相当于从i到j的总流量：[math]\displaystyle{ t_{i,j} }[/math]。

如果从i到j的所有路径中不存在着循环，如下图所示的树状结构：

我们很容易计算出从1到4的总流量为：

[math]\displaystyle{ t_{1,4}=50\times \frac{3}{5}\times \frac{2}{3}=20=T_1 (M+M^2)_{1,4} }[/math]

即1节点的总流量乘以经过一步概率转移从1到4的概率，再加上总流量乘以经过2步概率转移从1到4的概率。然而，当网络中存在着循环结构的时候，计算稍微麻烦，但是我们不妨将上面通过马尔科夫链计算的公式推广到有环的结构上，也就是对M的m次方求和一直进行下去，直到无穷大。例如，对于如图1所示的流网络，从1到4的总流量为：

[math]\displaystyle{ t_{1,4}=T_1 (\sum_{k=1}^{\infty}M^k)_{1,4} }[/math]

这里，T₁为从源到1的流量，等于80。我们知道M_1,4表示粒子直接从1到4的转移概率，这样T₁M_1,4就是粒子经过一步转移从1到4的流量。同样，[math]\displaystyle{ (M^2)_{1,4} }[/math]表示粒子经过两步转移从1到4的转移概率，那么[math]\displaystyle{ T_1(M^2)_{1,4} }[/math]就是粒子经过2步转移从1到4的流量。……，我们将经过1步转移、2步转移、……从1到4的粒子流总量加起来自然就是从1到4的粒子流量。

基础矩阵

然而，上式不方便计算，因为存在着矩阵的无穷级数和。注意到，M是一个行列式小于1的方阵，因此，我们可以用下列矩阵级数的求和公式：

[math]\displaystyle{ (\sum_{m=1}^{\infty}M^m)_{1,4}=\frac{M}{I-M} }[/math]

这里I是与M同样维度的单位阵。矩阵的无穷等比数列之和与数的无穷等比数列和的计算一模一样。我们定义：

[math]\displaystyle{ U=\frac{1}{I-M}=(I-M)^{-1} }[/math]

为基础矩阵（Fundamental Matrix），这沿袭了投入产出分析的用法。[math]\displaystyle{ (I-M)^{-1} }[/math]表示矩阵I-M的逆。在这个例子中，U可以写出来：

[math]\displaystyle{ \left( \begin{array}{ccccccc} 1 & 1 & \frac{3}{4} & \frac{7}{16} & \frac{1}{4} & \frac{1}{2} & 1 \\ 0 & 1 & \frac{3}{4} & \frac{7}{16} & \frac{1}{4} & \frac{1}{2} & 1 \\ 0 & 0 & \frac{42}{41} & \frac{7}{82} & \frac{14}{41} & \frac{28}{41} & 1 \\ 0 & 0 & \frac{12}{41} & \frac{42}{41} & \frac{4}{41} & \frac{8}{41} & 1 \\ 0 & 0 & \frac{3}{164} & \frac{21}{328} & \frac{165}{164} & \frac{21}{41} & 1 \\ 0 & 0 & \frac{3}{82} & \frac{21}{164} & \frac{1}{82} & \frac{42}{41} & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{array} \right) }[/math]

总流量

于是，将数字代入，从1到4的总流就是：

[math]\displaystyle{ t_{1,4}=T_1 (MU)_{1,4}=80\times \frac{1}{5}=16 }[/math]

下面，我们再来计算从2到4的总流量。按照上面的类似计算，我们可以得到：

[math]\displaystyle{ t_{2,4}=T_2 (MU)_{2,4}=60\times \frac{14}{41}=20.4878 }[/math]

然而，从图1中，我们看到从2到4的路径只有一条2-->4，所以从1到4的流量直观地看应该是20，而我们上面计算的结果却比20大。为什么呢？仔细分析会发现，原因出在，对于2号节点的总流量T₂包含着从2到4又循环回2的流量（从3到2的那部分流量）。而这部分流量又在U矩阵中被重复计算了。所以，要把它刨除掉。即只计算从0节点经过1流过来流量50，该流量就称为从源0到2的首达流。如果所有粒子都是无色的，而只有经过2号节点的粒子被染成红色，并且永不褪色，那么，从0到2的首达流就是：

[math]\displaystyle{ \phi_{0,2}=50+30\times \frac{2}{7}=\frac{410}{7} }[/math]

不难看出，无色粒子可能从1直接到2，也可能经过3到2。至于首达流的更一般的计算将放到下一节讨论。

[math]\displaystyle{ t_{2,4}=\phi_{0,2} (MU)_{2,4}=\frac{410}{7}\times \frac{14}{41}=20 }[/math]

这样的计算才可以得到正确的结果。因此，任意点i到j的总流量应该按下式计算：

[math]\displaystyle{ t_{ij}=\phi_{0,i} (MU)_{ij} }[/math]

根据上述公式，从i到i自己的总流量也就应该是：[math]\displaystyle{ t_{ii}=\phi_{0,i} (MU)_{ii} }[/math] 然而，这个式子并不完全。因为，它并不包括流[math]\displaystyle{ \phi_{0i} }[/math]自身。习惯上，我们定义从i到i的总流应该包括[math]\displaystyle{ \phi_{0i} }[/math]。 因此，写完全应该是：

[math]\displaystyle{ t_{ii}=\phi_{0,i} ((MU)_{ii}+1) }[/math]。

注意到，根据U的定义，[math]\displaystyle{ MU=U-I }[/math]，所以

[math]\displaystyle{ t_{ii}=\phi_{0,i} (U)_{ii} }[/math]。

最后，综合上述所有讨论，从i到j的总流量应该按照下式计算：

[math]\displaystyle{ t_{ij}=\phi_{0,i}(U)_{ij} }[/math]

即使对于i不等于j的情况上式也成立，因为U的非对角元与MU的非对角元相同。

首达流

上面的讨论中曾引出了首达流的概念，下面我们就来讨论一般地从i到j的首达流应该怎样计算。首先，我们定义从i到j的首达流为所有从i出发首次到达j的粒子总流量。还是假设所有的粒子经过i节点后被染成红色，并且同时假设粒子一旦经过j后，我们就把红色涂掉，那么每一时刻到达j的红色粒子总数就是从i到j的首达流。这样，所有经过i到达j的粒子就可以分成2类：'首达流和循环流，它们的总和就是总流动。

例如，对于图1来说，从1到4的总流动可以分为首达流和循环流：

[math]\displaystyle{ t_{1,4}=\phi_{1,4}+\psi_{1,4} }[/math]

我们已经知道t_1,4的计算了，如果我们会计算ψ_1,4的话，就很容易求出首达流了。所谓的从1到4的循环流一定是那些曾经至少已经经过一次4号节点的粒子所构成的流动。而我们知道，所有已经经过4号节点的粒子一定是那些首次到达了4号节点的粒子又重新从4出发，沿着所有可能路径再流回到4节点的粒子。这些粒子的路径就是如eq1式中所描述的，并且t不为0。同样的，根据流量与路径的关系，我们可以计算出：

[math]\displaystyle{ \psi_{1,4}=\phi_{1,4}(\sum_{m=1}^{\infty}M^m)_{4,4} }[/math]

也就是1到4的首达流再经过所有可能路径回到4的总流量。这样，我们就有：

[math]\displaystyle{ t_{1,4}=\phi_{1,4}+\psi_{1,4}=\phi_{1,4}+\phi_{1,4}(\sum_{m=1}^{\infty}M^m)_{4,4}=\phi_{1,4}(1+(\sum_{m=1}^{\infty}M^m)_{4,4})=\phi_{1,4}u_{4,4} }[/math]

其中，u_4,4表示(U)_4,4 这样就可以求出1到4的首达流：

[math]\displaystyle{ \phi_{1,4}=\frac{t_{1,4}}{u_{4,4}} }[/math]

我们定义，当i=j的时候：

[math]\displaystyle{ \phi_{i,i}=T_i }[/math]

而我们已经知道了t_1,4的表达式，因此1到4的首达流可以计算出，等于16*165/164=660/41。

上述方法适用于更一般的网络，因此，我们可以得到计算任一点i到j首达流的正确计算公式为：

[math]\displaystyle{ \phi_{i,j}=\frac{t_{i,j}}{u_{j,j}} }[/math]

小结

总结来看，我们这里定义了从i到j的两种流动：i到j的总流和i到j的首达流。如果将经过i节点的粒子染色，则总流就是每个时刻j节点接收到的红色粒子总数，而首达流就是所有那些第一次经过j节点的红色粒子总数。它们的计算方式如下：

总流：

[math]\displaystyle{ t_{i,j}=\phi_{0,i}u_{i,j}=\frac{t_{0,i}}{u_{i,i}}u_{i,j}=T_0\frac{u_{0,i}u_{i,j}}{u_{i,i}} }[/math]

首达流：

[math]\displaystyle{ \phi_{i,j}=\frac{t_{i,j}}{u_{j,j}}=T_0 \frac{u_{0,i}u_{i,j}}{u_{i,i}u_{j,j}} }[/math]

其中T₀为从源输入到整个网络的流量，也就是IS（参看流动网络）。根据这两个公式，可以计算任意网络任意节点对的总流和首达流。

我们看到，无论是总流还是首达流，都要用到U矩阵，而

[math]\displaystyle{ U=I+M+M^2+\cdot\cdot\cdot+M^{\infty} }[/math]

相当于包含了所有从i到j的一步转移概率、2步转移概率、3步转移概率……。但注意，u_ij并不是概率，因为u可能大于1。即便它不大于1，也要注意到M、M², ……矩阵每一行都是归一化的，所以把它们加起来就不再满足归一化条件（无论对于行还是列）。所以，u_ij并不表示从i到j的概率，而是表示i对j的总影响(Impact)。这是一个很含糊的概念，但目前也只能这样来解释U了。

而当u_ij乘以了一个流量φ_i，则所得的矩阵的每一个元素就有了明确的意义：表示从i到j的总流量。这个流量虽然从定义上看是一个瞬时t，i到j的粒子总数，但是根据计算式，它实际上是把t时刻的0步转移、t-1时刻的一步转移、t-2时刻的2步转移，……，所有时间步的粒子数加起来。由于流网络是平衡的，所以马尔科夫链是平稳的，所以我们可以把瞬时空间上的流量转变成所有历史的求和，这一点是最让人吃惊的事实。

平均时间（距离）

前面研究的是流量，下面我们研究与时间相关的主题。大量的粒子从网络中流动，我们假设单位时间粒子沿着有向连边跳转。这样，粒子需要m个时间步才能走完一条m长的路径。下面，我们将研究从网络上任意一点i到任意一点j的时间（距离）。

首先，由于网络中存在着环路，因此，从i流向j的粒子有可能经历无限长的时间步。为了避免这个问题，我们只考虑首达时间，也就是粒子从i出发第一次到达j的时间。但是，由于粒子可能沿不同的路径到达j，而这些路径的长度又不一样，所以没有确定的首达时间值。于是，我们只能用这些粒子的首达时间值来作为我们研究的目标，也就是平均首达时间。有趣的是，对于平衡的流网络来说，平均首达时间的计算存在着简洁而明确的表达式和计算方法。

任意点i到汇的平均时间

首先，我们来求从任意节点i（不包括汇）出发，到达汇N+1的粒子们的平均首达时间。由于到达N+1的粒子不再流向任何其它节点，因此这个平均首达时间其实就是粒子从i到达N+1的平均时间。下面，我们来计算从i到汇的平均时间。

我们不妨假设从任意节点i到达汇N+1的平均时间为l_i，那么，我们可以列出如下等式：

[math]\displaystyle{ l_i=1+\sum_{j=0}^{N}{m_{ij}l_j}, ~~~~\forall i\neq N+1 }[/math]

(eq2)

下面来解释这个式子。我们知道任何粒子要想跳到汇都有两种可能性：直接跳到N+1，或者先跳到任意一个其他的非汇节点j，然后再跳到汇N+1。对于第一种情况，它的路径长度确定性的是1，发生概率是[math]\displaystyle{ m_{i,N+1} }[/math]。而对于第二种情况，假设粒子跳到了j，再由j跳到N+1，那么我们知道粒子从j跳到N+1的平均首达时间是l_j，而粒子跳转到j的概率是m_i,j，因此粒子跳到任意的非汇节点，再跳到汇的平均路径长度就是[math]\displaystyle{ \sum_{j=0}^{N}{m_{ij}(1+l_j)} }[/math]，其中和号里的式子加1的原因是每条路经长度都要包含从i跳到j的这一步。于是，把这两种可能性合起来就有：

[math]\displaystyle{ l_i=m_{i,N+1}\cdot 1+\sum_{j=0}^{N}{m_{ij}(1+l_j)}=m_{i,N+1}+\sum_{j=0}^{N}m_{ij}+\sum_{j=0}^{N}m_{ij}l_j=1+\sum_{j=0}^{N}m_{ij}l_j, ~~~~\forall i\neq N+1 }[/math]

最后用到了马尔科夫链的归一化条件：[math]\displaystyle{ \sum_{j=0}^{N+1}m_{ij}=1 }[/math]。注意到eq2实际上是一个线性方程组，未知数有N+1个，方程个数也有N+1个，可以直接求解。事实上，如果我们定义：向量[math]\displaystyle{ L=(l_0,l_1,\cdot\cdot\cdot,l_N)^T }[/math]，则eq2可以写成矩阵的形式：

[math]\displaystyle{ L=\mu+M_{-(N+1)} L }[/math]

其中[math]\displaystyle{ \mu=(1,1,\cdot\cdot\cdot,1)^T }[/math]，1构成的N+1维列向量。[math]\displaystyle{ M_{-(N+1)} }[/math]表示原始马尔科夫矩阵去掉最后一行和最后一列。则，最终L可以解出：

[math]\displaystyle{ L=(I-M_{-(N+1)}^{-1})\cdot \mu=U_{-(N+1)}\cdot \mu=(\sum_{j=0}^{N}u_{ij})_{(N+1)\times 1} }[/math]

(eq3)

其中矩阵U_-(N+1)为U矩阵去掉最后一行和最后一列。最后一个等式是根据马尔科夫链删除节点的性质中的定理1得到的。这样L中的第i个元素就是i节点到汇的平均时间。

对于图1所示网络，我们可以计算出[math]\displaystyle{ L=(63/16, 47/16, 175/82, 66/41, 525/328, 197/164) }[/math]，其中从源到汇的平均路首达时间是63/16，也就是平均每个粒子经过3.9375步跳出整个网络。它刚好是[math]\displaystyle{ \sum_{j=0}^{N}u_{0j}=\frac{C_0}{T_0}+1=\frac{TST}{IS}+1 }[/math]。这后两个等式需要参考流网络的异速标度律.

从源到任意点i的平均首达时间

我们无法简单套用公式eq3来计算从源到任意点i的平均首达时间，这是因为如果将i节点去除，而非N+1，那么整个网络就有可能不连通，使得最后的求解出现问题。但是，我们仍然可以用类似的推理方法：我们假设0到任意i的平均首达时间为l_i。从源0到任意节点i的路径总可以分成两种情况，要么从0直接到i，要么从0到某个j再到i，于是假设0到j的距离为l_j，则第二种情况下的平均首达时间就是：[math]\displaystyle{ \sum_{j=1}^{N+1} m'_{ji}(l_j+1) }[/math]。注意，在这里，马尔科夫转移概率是逆向马尔科夫转移概率m'_ji而不再是m了，这是因为我们已经确定性地知道了粒子肯定到达i的情况下，它可能从j节点过来的概率。于是，在这种情况下，概率归一化的对象是给定在i观察到粒子的情况下，从所有可能的中间节点过来的概率。所以，如果使用m_ji的话，则不能满足概率归一的条件，所以只能用m'_ij。最后，再综合第一种情况，就有：

[math]\displaystyle{ l_i=m'_{0,i}+\sum_{j=1}^{N+1} m'_{ji}(l_j+1)=m'_{0,i}+\sum_{j=1}^{N+1} m'_{ji}+\sum_{j=1}^{N+1} m'_{ji}l_j=1+\sum_{j=1}^{N+1}m'_{ji}l_j,~~~~\forall i\neq 0 }[/math]

同样可以把它写作一个矩阵方程：

[math]\displaystyle{ L=\mu+M'_{-0}L }[/math]

于是，求解得：

[math]\displaystyle{ L=(I-M'_{-0})^{-1}\mu=U'_{-0}\mu }[/math]

注意，与eq3相比，我们把M换成了逆向马尔科夫矩阵M'。也就是说，求从源0到任意点i的首达时间问题相当于是把整个流网络倒转过来，把源当成汇，把汇当成源，然后求任意节点i到汇的首达时间问题。

对于图1所示网络，可以计算出： [math]\displaystyle{ L'=(1, 365/164, 193/82, 529/164, 285/82, 63/16) }[/math] 注意，第一个元素对应的是0到1的平均首达时间，……，以此类推。那么0到汇的平均首达时间是63/16，这与上面的计算结果是一致的。

任意节点i到任意节点j的平均首达时间

更一般地，我们应如何求出任意点i到j的平均首达时间呢？首先，我们需要明确这个平均首达时间到底是什么意思。还是假设所有流经i点的粒子都被染上了红色，并且粒子一经过j点就被抹掉红色。这样，j点接收的粒子就是首次到达j点的粒子。同时，假设每个粒子都有一个内置的秒表。它一旦被染成红色，秒表就开始计时，每个时间步走动一次。于是我们在j点读出粒子的计时，看它从i出发走了多少时间步过来。由于每个粒子的时间不一样，于是，我们就可以求出一个均值。这就是我们要求的平均首达时间。

如何计算这个平均首达时间呢？上面介绍的方法只适用于源或汇，而不能简单地拓展到非源或汇节点。下面，我们将介绍一种非常简单实用的方法，我称其为级数法。首先，根据平均首达时间的定义，我们知道从i到j的平均首达时间可以写为：

[math]\displaystyle{ l_{i,j}=\sum_{k=1}^{\infty}kp^{k}_{i,j} }[/math]

这里，p^k_i,j是粒子经过k步从i首达j的概率。这个概率怎么求呢？一种直观的想法就是用(M^k)_i,j来计算这个数值。但是，这是错误的，因为，p^k_i,j这个概率是指所有那些从i到j首达的粒子中，刚好经过了k步的粒子比例。而(M^k)_i,j表达的则是，在所有从i出发经过k步转移后的粒子中，那些刚好到达节点j的粒子比例。这两者显然不一样，因为从i出发的粒子经过k步后不见得都到达j。

让我们只考察到达j的那些粒子，每一时刻，这个流量是固定的，有一些粒子是从i出发首达的粒子，而其中有一部分刚好经历了k步。那么，从i到j的首达粒子就可以用首达流来计算。刚好经过了k步的首达粒子如何计算呢？我们知道Φ_0,i(M^k)_i,j是所有经过k步转移，从i到j的粒子数。但是这些粒子并不是首达粒子，有一些曾到过j又从j到j的粒子数。为了避免这种重复，我们不妨在M上做一些手脚：我们把M中j节点对应的行都设为0，即，定义：

[math]\displaystyle{ (M_{-j})_{p,q}=\left\{\begin{array}{ll} m_{p,q} & \mbox {if } p\neq j, \\ 0 & \mbox {if } p = j\end{array}\right. }[/math]

这样，概率p^k_i,j就可以这样计算：

[math]\displaystyle{ p^{k}_{i,j}=\frac{\phi_{0i}(M_{-j}^k)_{i,j}}{\phi_{i,j}} }[/math]

于是平均首达时间就变成了求下列无穷级数和：

[math]\displaystyle{ l_{i,j}=\sum_{k=0}^{\infty}k\frac{\phi_{0i}(M_{-j}^k)_{i,j}}{\phi_{i,j}}=\frac{\phi_{0i}(\sum_{k=1}^{\infty}kM_{-j}^k)_{ij}}{\phi_{i,j}} }[/math]

那么，下面的任务就是如何求解无穷级数[math]\displaystyle{ \sum_{k=1}^{\infty}k(M_{-j})^k }[/math]的和。为此，我们设此和为X，则在等式两边右乘以M_-j得到:

[math]\displaystyle{ XM_{-j}=\sum_{k=1}^{\infty}kM_{-j}^{k+1}=M_{-j}^2+2M_{-j}^3+\cdot\cdot\cdot }[/math]

与等式[math]\displaystyle{ X=M_{-j}+2M_{-j}^2+3M_{-j}^3+\cdot\cdot\cdot }[/math]相减得到：

[math]\displaystyle{ XM_{-j}-X=-(M_{-j}+M_{-j}^2+M_{-j}^3+\cdot\cdot\cdot)=-M_{-j}U_{-j} }[/math]

于是，解出：

[math]\displaystyle{ X=M_{-j}U_{-j}^2 }[/math]

其中，[math]\displaystyle{ U_{-j}=\sum_{k=0}^{\infty}M_{-j}^k=(I-M_{-j})^{-1} }[/math] 所以，i到j的平均首达时间就是：

[math]\displaystyle{ l_{i,j}=\frac{\phi_{0i}(M_{-j}U_{-j}^2)_{ij}}{\phi_{i,j}}=\frac{\phi_{0i}(M_{-j}U_{-j}^2)_{ij}}{\frac{\phi_{0i}u_{ij}}{u_{jj}}}=\frac{u_{jj}(M_{-j}U_{-j}^2)_{ij}}{u_{ij}} }[/math]

又根据马尔科夫链删除节点的性质，我们可以把上式简化为：

[math]\displaystyle{ l_{i,j}=\frac{u_{jj}(M_{-j}U_{-j}^2)_{ij}}{u_{ij}}=\frac{(MU^2)_{ij}}{u_{ij}}-\frac{(MU^2)_{jj}}{u_{jj}} }[/math]

于是我们得到了从任意节点i到j的平均首达时间计算公式。将这个公式套用到图1所示的网络上，我们可以得到从任意点i到j的平均首达时间，这就构成了平均首达时间矩阵：

[math]\displaystyle{ L= \left( \begin{array}{ccccccc} 0 & 1 & 88/41 & 373/164 & 7218/2255 & 557/164 & 63/16 \\ \infty & 0 & 47/41 & 209/164 & 4963/2255 & 393/164 & 47/16 \\ \infty & \infty & 0 & 9/4 & 58/55 & 5/4 & 175/82 \\ \infty & \infty & 1 & 0 & 113/55 & 9/4 & 66/41 \\ \infty & \infty & 3 & 2 & 0 & 1 & 525/328 \\ \infty & \infty & 2 & 1 & 168/55 & 0 & 197/164 \\ \infty & \infty & \infty & \infty & \infty & \infty & 0 \end{array} \right) }[/math]

我们看到，这与前面计算的从0到所有点，以及从所有点到汇的距离是一致的，这进一步验证了这个方法的正确性。

平均时间

根据上述级数法求解平均首达时间的启发，我们还可以求出从任意点i到达任意点j的平均时间而非平均首达时间。这里，平均时间按照如下方式定义：所有经过i的节点染成红色，并且永不退色，而且，经过i的粒子就开启一个计时器，则从j点处统计红色粒子，计算它们计时器所经历的时间。

有趣的是，按照上述方法，这个平均时间比平均首达时间更容易计算。还是按照平均时间的定义：

[math]\displaystyle{ k_{i,j}=\sum_{k=1}^{\infty}kp^{k}_{i,j} }[/math]

这里[math]\displaystyle{ p^{k}_{i,j} }[/math]的计算不用再考虑首达问题，于是

[math]\displaystyle{ p^{k}_{i,j}=\frac{\phi_{0i}(M^k)_{i,j}}{T'_{i,j}} }[/math]

这里，[math]\displaystyle{ T'_{i,j} }[/math]应为从i到j的总流。理应代入前面得到的关于[math]\displaystyle{ T_{i,j} }[/math]的表达式。这大体上正确，只有在i=j的时候不正确。因为在[math]\displaystyle{ T_{i,i} }[/math]中包含了从0到i的首达流[math]\displaystyle{ \phi_{0,i} }[/math]，而在平均时间的计算中，这部分流量不应包含，故而[math]\displaystyle{ T'_{i,i}=T_{i,i}-1 }[/math]。所以：

[math]\displaystyle{ T'_{i,j}=\phi_{0,i}(MU)_{ij} }[/math]

从而时间：

[math]\displaystyle{ p_{i,j}^k=\frac{(M^k)_{i,j}}{(MU)_{ij}} }[/math]

所以，平均时间的计算公式为：

[math]\displaystyle{ k_{i,j}=\frac{(MU^2)_{ij}}{(MU)_{ij}} }[/math]

根据这个公式，计算图1所示网络，得到彼此平均时间的矩阵为：

[math]\displaystyle{ K= \left( \begin{array}{ccccccc} \infty & 1 & 365/164 & 193/82 & 529/164 & 285/82 & 63/16 \\ \infty & \infty & 201/164 & 111/82 & 365/164 & 203/82 & 47/16 \\ \infty & \infty & 273/82 & 191/82 & 177/164 & 109/82 & 175/82 \\ \infty & \infty & 177/164 & 273/82 & 341/164 & 191/82 & 66/41 \\ \infty & \infty & 505/164 & 341/164 & 669/164 & 177/164 & 525/328 \\ \infty & \infty & 341/164 & 177/164 & 505/164 & 273/82 & 197/164 \\ \infty & \infty & \infty & \infty & \infty & \infty & \infty \end{array} \right) }[/math]

有趣地是，对角线元素表达的是节点i到i的平均回归时间，它可以用于节点i的中心度刻画。

[math]\displaystyle{ k_{i,i}=(\infty,\infty,273/82,273/82,669/164,273/82,\infty) }[/math]

该数值越小，说明这个节点在整个网络中的流量循环也就越快。

另外，我们还可以通过比较平均时间矩阵和平均首达时间矩阵，会发现任意两点之间的平均时间都比平均首达时间要大一些，通过做这两个矩阵的差，我们得到：

[math]\displaystyle{ K-L= \left( \begin{array}{ccccccc} & 0 & 13/164 & 13/164 & 223/9020 & 13/164 & 0 \\ & & 13/164 & 13/164 & 223/9020 & 13/164 & 0 \\ & & 273/82 & 13/164 & 223/9020 &13/164 & 0 \\ & & 13/164 & 273/82 & 223/9020 &13/164 & 0 \\ & & 13/164 & 13/164 & 669/164 &13/164 & 0 \\ & & 13/164 & 13/164 & 223/9020 &273/82 & 0 \\ & & & & & & \end{array} \right) }[/math]

观察这个矩阵，我们发现，所有有值的列除了对角线拥有相同的数值。这表明任何一点到达某一个点的平均时间与首达平均时间之差均为常数，并且出于同一环上的不同节点具有相同的数值，例如2,3,5。

参考文献

历史上，经济学家最早给出了U矩阵（也称为Leontief矩阵）的定义，请参照投入产出分析。后来生态学家Finn将这种分析方法引入到了生态流网络中用于分析能量流网络。Patten、Ulanowicz等人将这套方法系统化，提出了Environ理论。随后，这些人又将这套方法进一步扩展以讨论能量流的循环与存储问题以及能量流的节律等问题。

将流网络建模成马尔科夫链不仅仅可以讨论流量、时间等问题，还可以讨论如电导经过某点中转的平均首达时间等问题。

相关词条

流网络

编者推荐

集智课程

《网络科学》系列课程

复杂网络的数学建模与应用 | 网络科学第三期

网络科学已经跨过了二十个年头迎来了第三个十年。其发展历程大致如下:

第一个十年，根据实际的需要，网络科学重点研究了网络的基本模型及其性质，如相关性、鲁棒性、社团性等。除了已有的随机网络模型外，提出的最著名模型有小世界网络和无标度网络两种模型。我们网络科学第一期课程讲读巴拉巴西网络科学的书就涉及这些内容。

第二个十年，在网络传播和网络同步研究的基础上，网络科学的重点转入了网络动力学的全面研究，以及网络动力学与网络结构的相互关系的探索。在这十年，网络传播、博弈、同步和控制都得到了长足的进步。这些是我们网络科学第二期课程讲读汪小帆、李翔、陈关荣老师的书的内容。

第三个十年，已经过去了几年。从国内外网络科学研究的发展趋势来看，各种各样的更复杂的网络模型和结构以及高阶相互作用动力学引起了人们的极大兴趣。为了回应这种迫切需求， 我们网络科学第三期课程将围绕复杂网络的数学建模与应用进行多角度的介绍.

本期课程顺应国内外网络科学研究的发展趋势，从基础、应用、Plus进阶、代码实践模块，多角度介绍更复杂更贴近现实的网络模型和结构以及高阶相互作用动力学。