基于马尔科夫过程亚稳性的复杂网络簇结构分析和识别方 法 技术领域 本发明属于模式识别和数据挖掘领域， 尤其涉及社会网、 万维网和生物网络等复 杂网络的分析。
背景技术 现实世界中的诸多系统都以网络形式存在， 如社会系统中的人际关系网、 科学家 协作网和流行病传播网， 生态系统中的神经元网、 基因调控网和蛋白质交互网， 科技系统中 的电网、 因特网和万维网等。由于这些网络具有很高的复杂性， 因此被称为 “复杂网络” 。与 小世界性和无标度性并列， 复杂网络簇结构 (CNCS) 是复杂网络最普遍和最重要的拓扑结 构属性之一， 具有同簇节点相互连接紧密、 异簇节点相互连接稀疏的特点。 CNCS 识别方法旨 在揭示出复杂网络中真实存在的网络簇结构。
CNCS 识别方法对分析复杂网络的拓扑结构、 理解其功能、 发现其隐模式和预测其 行为都具有十分重要的意义， 具有广泛的应用前景， 目前已被应用于恐怖组织识别、 组织机 构管理等社会网络分析， 新陈代谢网络分析、 蛋白质交互网络分析和未知蛋白质功能预测、 基因调控网络分析和主控基因识别等生物网络分析， 万维网社区挖掘和基于主题词的万维 网文档聚类， 搜索引擎， 空间数据聚类和图像分割， 关系数据分析等众多领域。
目前已存在多种 CNCS 识别算法， 按照所采用的基本求解策略， 它们中的大多数可 归属为两大类 ： 基于优化的识别方法和启发式识别方法。 前者将 CNCS 识别问题转化为优化 问题， 通过最优化预定义的目标函数来计算复杂网络的簇结构， 后者基于预定义的启发式 规则设计启发式算法。谱方法和局部搜索方法是两类主要的基于优化的 CNCS 识别方法。
谱方法和局部搜索方法是两类主要的基于优化的 CNCS 识别方法。谱方法将网络 聚类问题转化为二次型优化问题， 通过计算特殊矩阵的特征向量来优化预定义的 “截” 函 数。谱方法具有严密的数学理论， 已发展成聚类的一种重要方法 ( 称为谱聚类法 )， 被广 泛应用于图分割和空间点聚类等领域。但是， 针对 CNCS 识别， 谱方法的主要不足是 ： 需要 借助先验知识定义递归终止条件， 即谱方法不具备自动识别网络簇总数的能力。1970 年的 Kernighan-Lin 算法、 2004 年的快速 Newman 算法和 2005 年的 Guimera-Amaral 算法是三个 典型的基于局部搜索优化技术的 CNCS 识别算法。这类算法都包含三个基本部分 ： 目标函 数、 候选解的搜索策略和最优解的搜索策略， 但在具体实现上各不相同。 其后提出的具有代 表性的基于优化的 CNCS 识别算法有 ： 2008 年提出的基于极大似然的 CNCS 识别算法、 2009 年提出的改进谱方法和 2010 年提出的基于蒙特卡洛方法的 CNCS 识别算法。
采用优化方法识别出的网络簇结构完全取决于优化目标， 因此 “有偏” 的目标函数 会导致 “有偏” 的解 ( 即得到的网络簇结构和真实存在的网络簇结构不符 )。 值得注意的是， 包括以上提到的快速 Newman 算法和 Guimera-Amaral 算法在内， 很多基于优化的 CNCS 识别 方法都以最大化 2002 年 Newman 提出的 Q 函数作为优化目标。然而研究发现， Q 函数是有偏 的， 并不能完全准确地刻画真实的网络簇结构。例如， 对于基准测试数据 Karate 社会网络
而言， 其真实的网络簇结构对应的 Q 值是局部极大值， 而非全局最大值。2004 年， Guimera 等研究发现， 对于某些随机网络， 由于受到扰动的影响， 明显不好的网络簇结构却对应相对 较高的 Q 值。2007 年， Fortunato 和 Barthelemy 系统地研究了 Q 函数对识别精度的影响， 指出 ： 基于优化 Q 函数的 CNCS 识别算法倾向于找到粗糙的而不是精细的网络簇结构。这意 味着， 在多数情况下这类算法不能识别出网络中真实存在的全部网络簇。
2002 年的 MFC(Maximum Flow Community) 算法、 2002 年的 Girvan-Newman(GN) 算 法、 2004 年的 Wu-Huberman(WH) 算法和 2005 年的 CPM(Clique Percolation Method) 算法 等是典型的启发式 CNCS 识别算法。其后提出的具有代表性的启发式算法还包括 ： 2008 年 提出的层次聚类算法、 2009 年提出的基于信息论的 CNCS 识别算法、 2010 年提出的基于拉普 拉斯动力性的多粒度 CNCS 识别算法等。 这类算法的共同特点是 ： 基于某些直观的假设来设 计算法采用的启发式信息， 对于大部分网络， 它们能够快速地找一个近似最优解， 但无法从 理论上严格保证它们对任何的输入网络都能找到令人满意的解。
综上所述， 尽管已存在多种方法， 但都具有各自的局限性， CNCS 识别问题还远远 没有被很好解决， 集中体现在以下 2 个方面 ：
第一， 从理论上我们还没有客观地认识清楚网络簇结构的本质含义。目前我们还 无法回答类似如下的基本问题 ： 网络簇结构是怎么形成的？它与网络的其它性质有什么必 然联系？它与网络自身的哪些内在属性有关？因此， 现阶段我们不得不通过观察有簇网络 所展示出的 “外在” 现象去理解网络簇概念， 进而借助 “主观” 定义的目标函数或启发式规则 去刻画和识别 CNCS。如前分析， 基于这些目标函数或启发式规则的算法常常会导致 “有偏” 的计算结果， 并且采用不同的目标函数或启发式规则常常会计算出不同的 CNCS。 因此， 一个 基本问题是 ： 从网络的 “内在” 属性出发， 我们能否给出一种 “客观” 的理论模型去解释、 刻 画和识别 CNCS。
第二， 现有的 CNCS 识别算法都具有各自的局限性， 不能同时满足无偏、 计算速度 快、 识别精度高、 无监督 ( 即不依赖先验知识、 对参数不敏感 ) 等基本要求。通过定性和定 量的分析、 比较现有的主要算法后发现， 识别精度高的算法往往具有很高时间复杂性 ( 高 于 O(n2))， 而快速的识别算法往往以牺牲精度为代价并且需要较多的参数和先验知识。另 外特别需要指出的是， 如何在没有任何先验信息的情况下识别出真实的网络簇总数仍是一 个未解决的难题。因此， 如何设计出快速、 高精度和无监督的 CNCS 识别方法是当前最期待 解决的问题之一。 发明内容 本发明的目的在于揭示复杂网络分簇现象的本质， 并提供一种用 于定量分析和 快速识别复杂网络簇结构的方法。
为实现上述目的， 本发明提供了一种基于马尔科夫过程亚稳性的复杂网络簇结构 分析和识别方法， 其特征在于包括如下步骤 ：
构造给定复杂网络上的一个 Markov 过程 ；
计算该过程的一步转移概率矩阵， 并计算该矩阵的特征值 ；
通过分析特征值计算出网络簇个数 ；
计算 Markov 过程的第一个亚稳态 ；
根据第一个亚稳态识别出网络的全部网络簇及其层次结构。附图说明 图 1 所示的流程图给出了基于 Markov 过程亚稳性分析和识别复杂网络簇结构的 理论框架 NAP ；
图 2 所示的流程图给出了以上理论框架的一种快速实现方法 fast_NAP， 该方法能 够在近似线性的时间内识别出网络中真实存在的全部网络簇及其层次结构。该方法免参 数， 且不需要网络先验知识 ；
图 3- 图 7 给出了采用 NAP 和 fast_NAP 方法分析不同网络的结果。
具体实施方式
下面将对本发明进行详细说明。
参照图 1， NAP 方法的流程开始于步骤 101。
步骤 102 给出了构造网络上 Markov 过程的方法， 具体如下 ：
假设网络中存在一个 Agent， 该 Agent 能够沿着网络连接从一个网络节点随机的 移动到其它网络节点。用 P{Xt ＝ i， 1 ≤ i ≤ n} 表示 Agent 经过 t 时间到达网络节点 i 的 概率， X ＝ {Xt， t ≥ 0} 表示 Agent 在不同时刻位置构成的随机序列。由于 Agent 在 t 时刻 的位置唯一地由其在 t-1 时刻的位置决定而和 t-1 时刻之前的位置都无关， 即满足 ：
P(Xt ＝ it|X0 ＝ i0， X1 ＝ i1，…， Xt-1 ＝ it-1} ＝ P(Xt ＝ it|Xt-1 ＝ it-1}
因此随机序列 X 满足马尔科夫性， 是一个离散、 时齐的 Markov 过程。
步骤 103 给出了构造网络上 Markov 过程转移概率矩阵 P 的方法， 具体如下 ： -1
P＝D A
其中矩阵 A ＝ (aij)n×n 表示网络的邻接矩阵， D ＝ diag(d1，… dn) 表示网络的度 矩阵， 其中 di ＝∑ jaij.
步骤 104 给出了计算 I-P 矩阵特征值的方法， 具体采用幂法， 迭代地计算出矩阵 I-P 的全部特征值。
步骤 105 给出了根据以上特征值计算出网络簇个数 K 的方法， 具体如下 ：
其中， CQK ＝ λK/λK+1， Λ ＝ (λ1，…， λn) 表示矩阵 I-P 的特征值。 步骤 106 给出了计算出 Markov 过程第一个亚稳态 S1 的方法， 具体如下 ：其中 P 表示步骤 103 计算出来的 Markov 过程一步转移概率矩阵。λK+1 表示矩阵 I-P 的第 K+1 小特征值。
步骤 107 给出了从 S1 识别出 K 个网络簇及其层次结构的一般方法， 其中一个具体 方法如下 ：
将矩阵 S1 中的每一行看作一个 n 维的向量， 于是共有 n 个 n- 维向量， 每个向量对 应网络中的一个节点 ；
采用基于相似度的向量聚类方法 ( 如 K- 均值 ) 将以上 n 个 n- 维向量聚成 K 个类，
对应网络的 K 个簇。
参照图 2， fast NAP 方法的流程开始于步骤 201 ；
步骤 202 给出了选择网络中最稳定的状态的方法， 具体如下 ：
其中 di 表示节点 i 的度。 步骤 203 给出了计算 c 的不变序分布的方法， 具体如下 ： 迭代计算所有状态到达状态 c 的 t- 步转移概率 ：按照概率值将所有状态排序， 得到 t 时刻状的态索引序列 St ； 为c的重复以上两步直到状态索引序列不变为止， 即 St ＝ St-1 ； 状态索引序列不变时刻 T 对应的状态转移概率分布不变序分布 OTD.
步骤 204 给出根据 c 的不变序分布计算出网络的一个最优二分的方法， 具体如下： 根据计算出的不变序分布 OTD， 将网络节点按照 OTD 各分量取值由小到大排序 ；
按照以上节点次序将网络 N 的邻接矩阵 A 转换为新的邻接矩阵 B ；
令 I{f} 表 示 试 性 函 数， 如 果 f 为 真 则 I{f} 取 值 1， 否 则 取 值 0. 令 向 量 Xx(i) ＝ 1·I{i ≤ x}+(-1)·I{i ＞ }， 1 ≤ i ≤ n， 表示网络划分 (N1， N2)， x(1 ≤ x ≤ n) 表示网络节点序 的分割点， 通过 x 可将节点序分成两个部分， 进而将整个网络 N 分成两个子网络 N1 和 N2. 最 * 优网络划分 x 满足 ：
其中向量 Yx ＝ (1-x)(1+Xx)-x(1-Xx)， D 表示步骤 103 计算出的网络度矩阵， QB ＝I-D B.
-1根据上式， 可以依次计算出 x ＝ 1， x ＝ 2， ...， x ＝ n 对应的 n 个值， 其中最小值对应的 x 值即为待求的网络最优划分 x*， 进而将网络 N 划分为 N1 和 N2 两个子网 络。
步骤 205 给出停止条件的判断方法， 具体如下 ：
EP(N1， N2) ≥ 0.5 ∧ EP(N2， N1) ≥ 0.5
其中，表示从 N1 到 N2 的逃逸概率， (N1， N2) 表示网络 N 的一个二分， 满足 N1 ∪ N2 ＝ N 且 N1 ∩ N2 ＝ Φ， 其中 Φ 表示空集。P 是由步骤 103 计算出的 Markov 过程转移概率矩阵。
步骤 207 给出将网络划分为两个子网络方法， 具体如下 ： *
根据步骤 204 计算出来的最优分割点 x 为向量 Xx 赋值 ：Xx 分量值为 1 对应的节点及其相互间的连接构成第一个子网络， Xx 分量值为 -1 对应的节点及其相互间的连接构成第二个子网络。
步骤 208 给出递归处理第一个子网络的方法， 具体如下 ： 将 207 步骤得到的第一个 子网络作为输入从步骤 202 开始执行。
步骤 209 给出递归处理第二个子网络的方法， 具体如下 ： 将 207 步骤得到的第二个 子网络作为输入从步骤 202 开始执行。
以下采用不同的网络对 NAP 和 fast_NAP 方法进行了测试与评价， 进一步说明本发 明方法的原理和效果， 具体如下 ：
例 1 采用 NAP 方法分析网络中真实存在的网络簇个数
图 3(a) 表示一幅包括 4 个字符图像。图 3(b) 表示该图像对应的网络， 网络建模 方法采用全连接法， 网络连接的权值采用高斯相似度公式计算。 该网络中， 每个字符构成一 个自然的网络簇， 因此该网络的真实网络簇个数为 4。图 3(c) 表示采用 NAP 计算出来的各 个 CQK 值， 其中 CQ4 最小， 因此 NAP 计算出的网络簇总数为 4， 和真实的网络簇总数一致。
例 2 图 4 给出了实现 fast_NAP 方法的软件界面。采用 fast_NAP 方法， 该软件能 够识别网络中的全部网络簇及其层次结构， 并采用邻接矩阵和层次树相结合的方式可视化 复杂网络簇及其层次结构。 通过重新安排原始邻接矩阵的行和列， 将同簇节点排列在一起， 可以得到能够清晰表示网络簇结构的转换邻接矩阵。如果网络有分明的簇结构， 其对应的 转换矩阵应该是一个近似对角矩阵， 主对角线的每一个块子矩阵 恰好对应一个网络簇。 分 布在主对角线区域的非零元素 ( 对应簇内边 ) 远远多于散落在主对角线区域之外的非零元 素 ( 对应簇间边 )。簇间的层次关系用树形表示。
例 3 采用 fast_NAP 方法分析美国大学足球盟网络
图 5(a) 给出了 2000 赛季美国大学足球联盟网络。该网络包含 115 个节点和 613 条边。网络中每个节点表示一支大学足球队， 每条边表示两个球队间进行的一场比赛。根 据地理位置， 全部球队组织成 12 个联盟。根据比赛规则， 联盟内的比赛远远多于联盟间的 比赛。因此， 按照比赛的关系， 12 个联盟对应 12 个网络簇。
图 5(b) 给出了 fast_NAP 方法的计算结果， 一个对角化的邻接矩阵和一个网络簇 层次结构树。分析比较后发现， fast_NAP 算法得到的 12 个网络簇与 12 个实际的足球联盟 基本吻合， 只有 6 支隶属于 3 个相对独立的联盟的球队被错分。错分的原因在于， 这些球队 与联盟外球队的比赛过多。
例 4 采用 fast_NAP 方法分析海豚网络
图 6(a) 给出了海豚社会网络。该网络包含 62 个节点， 160 条边。每个节点表示一 只海豚， 每条边表示两只海豚之间的社会关系。这群海豚被跟踪观察 7 年。由于某种原因， 该群海豚最终分裂为 2 个子群， 如图 6(a) 所示。
图 6(b) 给出了 fast_NAP 方法的计算结果， 一个对角化的邻接矩阵和一个网络簇 层次结构树。分析比较后发现， fast_NAP 方法得到的 2 个最大网络簇准确预测了海豚社会 的实际分裂情况。 此外，fast_NAP 方法还分析出了分裂后的 2 个海豚群可能进一步发生的 分裂情况。
例 5 采用 fast_NAP 方法进行图像模式识别
本节借助图像模式识别问题测试 fast_NAP 方法的有效性。当图像被建模成网络 后， 图像中主要元素的识别就转换为相应的网络簇识别。图 7 的第一列给出了 3 幅不同图 像， 它们对应的网络按照以下方法生成 ：
(1) 每个像素看作是一个网络节点 ；
(2) 对于每一对像素 i(xi， yi) 和 j(xj， yj)， 根据高斯相似度公式计算出连接 (i， j) 的权值 ；
(3) 去除权值低于平均值的连接。
图 7 的第二列给出了按照以上方法得到的不同图像对应的网络。
图 7 的第三列给出了应用 fast_NAP 方法的识别结果， 其中不同的图像元素用不同 的颜色区分。如图所示， 所有图像元素都被正确识别。特别指出的是， 在第三幅图像中， 字 符 “A” 和 “B ” 能够从噪声中被正确识别。