基于图形处理器的大规模静态网络划分方法

摘要

本发明公开一种基于图形处理器的大规模静态网络划分方法，主要是针对现有网络划分方法不适于大规模网络的划分，且划分方法效率低而设计。本发明通过定义一个网络划分结果矩阵S，将网络划分中分别求各子模块的特征向量和特征值转换成求整个网络的特征向量和特征值，并根据图形处理器的并行计算体系结构，将网络划分中的特征值和特征向量的计算分解为几个基本运算，将原本庞大的稠密矩阵计算问题简化为稀疏矩阵和向量之间的运算。本发明通过基本运算间的合理组合节省了内存空间，并且减少了数据传输的消耗，有效地提高了网络划分的效率，将本不易实现的超大规模网络划分成为可能。

说明书

技术领域

本发明涉及计算机以及电子信息技术领域，尤其涉及大规模静态网络划分的方法。

背景技术

网络模型在近年来的统计物理和应用数学研究领域得到了极大的重视，而网络这一概念也被成功的应用到了其他领域，如互联网，社交网络以及生态网络等。网络研究包括网络各种属性的计算，包括小世界属性，重尾度分布和群落结构；其中，网络群落结构的研究对于网络研究具有重要的意义，比如对于社交网络群落结构的研究可以帮助发现群体的共同爱好和共同特点。网络群落结构研究可以帮助理解网络的结构功能以及拓扑关系。

研究网络群落特性的主要方法是网络划分。网络划分采用的是图论中的划分方法，图的划分在图论研究中占有主要的地位。目前比较成熟的划分方法包括最小切割法，层次聚类法，随机游动法和格文纽曼法。还有一种是基于特征值计算的划分方法，它能够得到较为准确的划分结果，并且具有充分的理论支持，下面是这种方法的详细内容：

这种网络划分方法针对无权无向网络。所谓无权无向网络，即网络中连接各节点的边没有方向和大小。无权无向网络可以由邻接矩阵表示，所谓邻接矩阵，即矩阵的元素表示节点的边的有无，有边则值为1，否则为0。如下式所示，A为网络的邻接矩阵，A_ij为邻接矩阵第i行第j列的元素。

为实现网络的划分，需定义网络模块化系数，该系数反映出网络划分的程度。当网络是随机网络时，模块化系数为零；当网络的划分越接近其本身具有的集群属性时，模块化系数越大，所以网络划分的目的在于最大化模块化系数。

模块化系数定义为网络每两点间实际连接数与连接数期望之差的和，所谓连接数期望即有网络各点属性计算出的两点间的预期边数。模块化系数Q定义如下：

$Q=\frac{1}{4m}\underset{\mathrm{ij}}{\Sigma }[A_{\mathrm{ij}}-P_{\mathrm{ij}}](s_i{s}_j+1)$

其中，A_ij为节点i与节点j间的连接数(0或1)，P_ij为两点间连接数期望，m为网络总边数。S_i用来决定网络的划分，当节点i与节点j属于同一个模块时s_is_j＝1，否则s_is_j＝-1。网络划分的目的即是寻找合适的S_i使Q值最大。连接数期望可以表示为：

$P_{\mathrm{ij}}=\frac{k_i{k}_j}{2m}$

其中，k_i为节点i的度，即与节点i连接的边数，m为网络总边数。

这一寻找过程可以转化为矩阵特征值计算的过程，下面以将网络划分成两部分为例说明如何通过求矩阵特征值的方式求出网络的划分方式。

首先，定义模块化矩阵B的元素B_ij：

B_ij＝A_ij-P_ij

当S_i对应为矩阵B的最大正特征值的特征向量时，Q取最大值。由于只将矩阵划分成两部分，所以S_i只可取+1或-1，而S_i具体取正还是取负要根据求出的特征向量各元素的正负来确定。综上所述，当将网络划分为两部分时，只需求出B矩阵最大特征值对应的特征向量即可。

将网络划分成多个模块，就是不断地将已划分的子块继续分成两块，由于不能忽略子块与整个网络其他部分的联系，不能完全将子块看作一个单独的矩阵。因此需要进行如下的处理：

设有一子块G，共有n_g个节点，对应一个n_g×n_g的模块化矩阵B^(G)：

$B_{\mathrm{ij}}^{\left(G\right)}=B_{\mathrm{ij}}-{\delta }_{\mathrm{ij}}\underset{l\in G}{\Sigma }{B}_{\mathrm{il}}$

其中，l为子块G的索引值，表示属于子块G的点的下标；而i和j也对应为子块G元素的索引。

如上所述，整个划分过程是自顶向下的，即先将整个网络划分为两个模块，再将每个模块继续划分为两个模块，以此类推。

这种方法能够提供较为准确的结果，并且具有严格的理论证明保证解为最优，可是方法的实现较为复杂，无法承受大规模网络的计算。方法的复杂性主要表现在如下两个方面：第一，涉及到模块化矩阵B^(G)的特征值的计算，而当网络维度很大时，矩阵特征值的计算量是巨大的，例如使用雅可比迭代算法时，每次迭代的时间复杂度为O(n³)；第二，在将网络不断二分的过程中，子块的大小会发生变化，并且无法保证子块下标的连续性，因此不利于算法的并行化。

发明内容

针对上述问题，本发明提供一种能够有效提高网络划分的效率，实现超大规模的网络划分的基于图形处理器的大规模静态网络划分方法。

为达到上述目的，本发明基于图形处理器的大规模静态网络划分方法，包括以下几个步骤：

(1)定义待划分网络的邻接矩阵A，各节点度向量K，网络划分结果矩阵S和模块化矩阵B’；

(2)中央处理器将已定义的各矩阵和向量传向图形处理器；

(3)图形处理器接收中央处理器传输的数据，计算出模块化矩阵B’的特征向量和特征值，并将其传输到中央处理器中；

(4)中央处理器接收上述特征向量和特征值，对特征向量进行二值化处理，得出两个划分向量，并将其分别赋值给网络划分结果矩阵S的相应行向量；然后，判断特征值是否大于零，是，继续重复步骤1～4；不是，输出网络划分结果矩阵S。

进一步地，所述的网络划分结果矩阵S为一n×m矩阵，

其中，n是待划分网络的节点数，m是待划分网络划分完成后的总模块数。

进一步地，步骤(1)具体包括以下步骤：

1.1、定义待划分网络的邻接矩阵A，各节点度向量K和网络划分结果矩阵S；

1.2、定义模块化矩阵B’：

B′＝B-diag{d₁，...，d_n}

B＝A-P

$P=\frac{k·k^T}{2m}$

$d_i=\underset{l}{\Sigma }{A}_{\mathrm{il}}-\frac{k_i}{2m}\underset{l}{\Sigma }{k}_l=A_i·s^k-\frac{k_i}{2m}{k}^T·s^k$

得，$B^{\prime }=A-\frac{k·k^T}{2m}-\mathrm{diag}\{d_1,...,d_n\}$

其中，P为边数期望矩阵，A_i为邻接矩阵A的第i行，S^k为网络划分结果矩阵S的第k行向量。

进一步地，步骤(3)还包括：定义一向量u_n，将模块化矩阵B’与向量u_n相乘，基于乘幂法计算出模块化矩阵B’的特征向量和特征值。

进一步地，步骤(3)具体实现步骤如下：

3.1、定义一n维向量u_n；

3.2、将模块化矩阵B’与向量u_n相乘，即：

$B^{\prime }·u_n=A·u_n-\frac{k^T·u_n}{2m}k-\mathrm{diag}(d_1,L,d_N)·u_n$

3.3、基于乘幂法计算，即：

$u_{n+1}=\frac{B^{\prime }·u_n}{\left|\right|B^{\prime }·u_n\left|\right|}$

当n→∞时，X≈u_n；β≈||u_n||

式中，X为模块化矩阵B’的主特征值对应的特征向量，β为模块化矩阵B’的主特征值。

特别地，步骤3.3还包括：通过至少两次计算后，依据u_n中各元素的正负变化确定主特征值的正负；若为负，将模块化矩阵B’与以主特征值为对角元的对角阵相加，再进行一次乘幂法求出主特征值。

进一步地，所述的两个划分向量，其中一个划分向量为特征向量中大于零的元素均为1，小于等于零的元素均为0的向量；另一个划分向量为特征向量中小于零的元素均为1，大于等于零的元素均为0的向量。

本发明具有以下几点有益的效果：

1、本发明通过定义一网络划分结果矩阵S，将原先求某一子块的模块化矩阵B^(G)的特征向量转化为求整个n维网络特征向量的问题。其次，本发明保证了内存的连续，在图形处理器中的并行计算性能得到了改善。

2、本发明将B’与向量u_n相乘分解为稀疏矩阵与向量乘积和向量与向量乘积，再将整个乘法分解为各个基本运算，充分地利用了图形处理器的并行计算处理能力，有效的节省了内存的使用空间，提高了求解的速度。

附图说明

图1为本发明基于图形处理器的大规模静态网络划分方法的原理图；

图2为本发明所述图形处理器内部计算原理框图；

图3为本发明基于图形处理器的大规模静态网络划分方法的计算程序流程图。

具体实施方式

下面结合说明书附图对本发明的具体实施方式做详细描述。

如图1所示，本发明基于图形处理器的大规模静态网络划分方法的原理图。如图中所示，中央处理器执行数据的调度和网络的划分，图形处理器执行特征值和特征向量的计算。整个网络划分过程的实现由两个数据队列实现，其一是待划分模块向量S^k队列，另一个是网络划分结果矩阵S各行向量队列；其中，网络划分结果矩阵S矩阵是一个n×m的矩阵，n是网络的节点数，m是总模块数。S的每一行均为一个n维的向量，表示一个模块，如下式：

其中，k表示模块序号，即S矩阵的某一行。

初始时，待划分模块向量S^k队列中只有一个向量S¹＝(1，1，...1)，即网络中的各节点都在一个模块中。中央处理器调用已定义待划分网络的邻接矩阵A，各节点度向量K，以及向量S¹＝(1，1，...1)；输入到图形处理器中。图形处理器求出模块化矩阵B’的特征值β和特征向量X，并将求出的特征值β和特征向量X传输回中央处理器。

中央处理器首先对特征向量X进行二值化处理，处理后得出两个划分向量X₁和X₂，其中，划分向量X₁为特征向量X中大于零的元素均为1，小于等于零的元素均为0的向量；划分向量X₂为特征向量X中小于零的元素均为1，大于等于零的元素均为0的向量。两个划分向量X₁和X₂即为输入模块向量S¹划分后的两个子模块。然后，将这两个子模块分别赋值给相应行向量，如S¹＝X₁，S²＝X₂；最后，判断计算结果特征值β，若特征值β＞0，表示划分后的两个子模块仍可再分，将向量S¹和S²存入待划分模块向量S^k队列中，再分别对模块S¹和S²按上述划分步骤进行进一步划分；若特征值β≤0，表示划分后的两个子模块S¹和S²已无法再分，输出网络划分结果矩阵S＝(S¹，S²)^T。当输入队列为空时，算法结束。

其中，模块化矩阵B’的定义如下：

B′＝B-diag{d₁，...，d_n}

$d_i=\underset{l}{\Sigma }{A}_{\mathrm{il}}-\frac{k_i}{2m}\underset{l}{\Sigma }{k}_l$

B＝A-P

$P=\frac{k·k^T}{2m}$

得，$B^{\prime }=A-\frac{k·k^T}{2m}-\mathrm{diag}\{d_1,...,d_n\}$

其中，A是网络邻接矩阵，k是各节点的度向量，P为边数期望矩阵，l是S^k为1位置的下标，d_i还可以进一步表示为矩阵A_i与S^k的乘积，即：

$d_i=A_i·s^k-\frac{k_i}{2m}{k}^T·s^k$

其中，A_i为邻接矩阵A的第i行，是一个n维向量。

B’矩阵的定义是由B和B^(G)的定义演化过来的，目的是将原本部分矩阵的划分转化为完整矩阵的划分，并且可以将乘幂法中原本稠密矩阵的运算转化成稀疏阵和向量间运算。B和B^(G)的定义在背景技术中有详细描述。

如图2所示，本发明所述图形处理器内部计算原理框图。图形处理器中的程序主要负责计算模块化矩阵B’的特征值β和特征向量X，主要功能可以分为预处理部分1和特征值特征向量求解部分2。预处理部分1用于计算d_i；特征值特征向量求解部分2通过乘幂法的修改形式计算模块化矩阵B’的最大正特征值和特征向量。特征值和特征向量计算的核心步骤是迭代计算B’·u_n，具体内容详述如下。

特征值和特征向量求解的过程中最困难也是最突出的问题就是模块化矩阵B’是一稠密阵，而稠密阵的特征值和特征向量的计算是十分复杂的。本发明为了提高特征向量和特征值的计算效率，使用图形处理器对其进行了改进。

首先，将模块化矩阵B’化简，即将B’与向量u_n相乘分解为稀疏矩阵与向量乘积，以及向量与向量乘积；然后，将整个乘法分解为各个基本运算；最后，采用乘幂法计算矩阵的特征值和特征向量。

其计算过程如下：

1、定义一n维向量u_n；

2、将B’与向量u_n相乘，即：

$B^{\prime }·u_n=A·u_n-\frac{k^T·u_n}{2m}k-\mathrm{diag}(d_1,L,d_N)·u_n$

3、基于乘幂法计算，即：

$u_{n+1}=\frac{B^{\prime }·u_n}{\left|\right|B^{\prime }·u_n\left|\right|}$

当n→∞时，X≈u_n；β≈||u_n||

式中，X为模块化矩阵B’的绝对值最大特征值对应的特征向量，β为模块化矩阵B’的主特征值，即绝对值最大特征值。图2中的这一步是乘幂法步骤中的一步，将结果矩阵归一化，以免出现当U_n+1过于庞大可能出现的数值越界情况，相当于

经过两次计算后，u_n中各元素的正负变化可以确定主特征的正负，如果主特征值为负，需将模块化矩阵B’与以主特征值为对角元的对角阵相加，此时矩阵全部特征值为正，再进行一次乘幂法求出主特征值，即为模块化矩阵B’最大正特征值与主特征值的和，减去主特征值后即可得到最大正特征值。

从上式可以看出，原本稠密矩阵与向量的乘法转化成三个计算的组合，其中A·u_n为稀疏矩阵与向量的乘法，为向量间乘法，diag(d₁，L，d_N)·U_n为向量数乘。

通过将原本稠密的矩阵向量乘转化为稀疏矩阵向量乘和向量乘法后使整个运算对内存的消耗大大减小，只需要存储稀疏阵和向量即可。另外，在分解计算后，每一步的划分只需更新diag(d₁，L，d_N)这部分对角矩阵，其余部分都不需改变，这样也大大减少了数据传输的消耗。

如图3所示，本发明基于图形处理器的大规模静态网络划分方法的计算程序流程图。图中各算式后面的标号表示基本运算的种类。

以下为各基本运算的具体说明：

(1)、稀疏矩阵与向量的乘法

A为稀疏矩阵，a为向量，乘法结果为向量b，有：

b＝A·a

(2)、向量间乘法

向量a与向量b相乘，结果为实数c：

c＝b^T·a

(3)、向量加法与数乘

向量a与向量b线性组合，结果为向量c：

c＝αa+βb

(4)、向量直乘

向量a与向量b直乘，结果为向量c：

$c=a\otimes b\Rightarrow c_i=a_i{b}_i,$i＝0...n-1

(5)、向量r范数

向量a的r范数为||a||_r：${\left|\right|a\left|\right|}_r=\sqrt[r]{\underset{i=0}{\overset{n-1}{\Sigma }}{a}_i^r}$

每种运算都能在图形处理器上得到加速，整体运算过程通过几种基本运算的组合也得到了简化，最终实现了整个算法性能的提高。

下面分别说明各基本运算在图形处理器上的实现方法。

图形处理器上各个基本算法是基于一种分合思想而实现的。在图形处理器中，“分”这一步是并行进行的，每一个线程计算两向量对应下标元素的乘法。“合”这一步是以二进制计数的方式在对应线程中将每一线程计算出的结果进行加法计算。

1.稀疏矩阵与稠密向量的乘法

稀疏矩阵与稠密向量的乘法可以分解成多个稀疏向量和稠密向量的乘积。其中，“分”运算为：图形处理器的每一个计算线程稀疏向量与稠密向量对应下标的元素的乘积；“合”运算为：在对应线程中将每一个线程计算出的乘积进行求和。

在本发明中，稀疏矩阵为邻接矩阵，所有值皆为1；且是以压缩稀疏行的形式存储的。于是，整个运算即转化为一个收集的过程，即根据稀疏矩阵的索引下标将稠密向量的对应元素进行相加运算。

所谓压缩稀疏行存储是把稀疏矩阵以稀疏行向量的形式存储。

2.向量间乘法

由向量乘法的数学表达式可以将向量乘法划分为乘法和加法两个部分，乘法可以看作“分”的运算，而加法可以看作“合”的运算。“分”的运算是将两个向量下标相同的元素相乘；“合”的运算是将相乘的结果相加。

“分”运算：π_i＝a_ib_i，i＝0...n-1

“合”运算：$c=\underset{i=0}{\overset{n-1}{\Sigma }}{\pi }_i$

3.向量间加法与数乘

在向量加法与数乘的运算中，“分”运算是两向量下标对应元素之间的乘法和加法，如下式所示：

“分”运算：c_i＝αa_i+βb_i，i＝0...n-1

由于输出为向量，没有“合”运算这一步，c_i即为输出向量的元素。

4.向量直乘

在向量加法与数乘的运算中，“分”运算与向量乘法相同，如下式所示：

“分”运算：c_i＝α_ib_i，i＝0...n-1

由于输出为向量，没有“合”运算这一步，c_i即为输出向量的元素。

5.向量r范数

在向量r范数的运算中，“分”运算为向量元素的r次幂运算，“合”运算为各“分”运算结果之和，如下式所示：

“分”运算：${\pi }_i=a_i^r,$i＝0...n-1

“合”运算：${\left|\right|a\left|\right|}_r=\sqrt[r]{\underset{i=0}{\overset{n-1}{\Sigma }}{\pi }_i}$

按照以上所介绍的各个基本运算的实现方法，本发明大大地提高了大规模网络划分的运算效率。下表为图形处理器实现的网络划分方法与中央处理器实现的加速结果：

  中央处理器耗时  图形处理器耗时  加速比 100次划分  74496.97  928.73  80

 最快划分  2515.55  31.68  79 最慢划分  3.78  0.29  13

本发明所述基于图形处理器的大规模静态网络划分方法的理论基础和具体实现步骤如下：

首先，数据结构的准备。

由于这种网络划分采用求特征值的方法，实际的划分向量S取值并不影响方法的实现。本发明以S矩阵来表示网络的划分结果。S矩阵是一个n×m的矩阵，n是网络的节点数，m是总模块数。S的每一行均为一个n维的向量，表示一个模块，如下式：

其中，k表示模块序号，k＝1，2...，即S矩阵的某一行。

这样便可以利用划分向量S^k，将求解n_g×n_g维B^(G)的特征向量的问题转换成求整个n维网络特征向量的问题。

然后，定义模块化矩阵B’：

B′＝B-diag{d₁，...，d_n}

$d_i=\underset{l}{\Sigma }{A}_{\mathrm{il}}-\frac{k_i}{2m}\underset{l}{\Sigma }{k}_l$

B＝A-P

$P=\frac{k·k^T}{2m}$

得，$B^{\prime }=A-\frac{k·k^T}{2m}-\mathrm{diag}\{d_1,...,d_n\}$

其中，A是网络邻接矩阵，k是各节点的度向量，P为边数期望矩阵，l是S^k为1位置的下标，d_i还可以进一步表示为矩阵A_i与S^k的乘积，即：

$d_i=A_i·s^k-\frac{k_i}{2m}{k}^T·s^k$

其中，A_i为邻接矩阵A的第i行，是一个n维向量。

依据上述算式，即可计算出模块化矩阵B’的特征值β和特征向量X。虽然根据B’计算得到的特征向量虽然是n维的，但是在不属于该模块的点的位置处的元素均为0，所以与计算B^(G)得到的特征向量是等价的。采用本发明的方法，虽然每次特征值计算中矩阵维数不会减少，但是由于保证了内存的连续，在并行计算中性能得到了有效地改善。

为了进一步提高提高特征向量和特征值的计算效率，可以将B’与向量U_n相乘分解为稀疏矩阵与向量乘积，以及向量与向量乘积；然后，将整个乘法分解为各个基本运算；最后，采用乘幂法计算矩阵的特征值和特征向量。

其计算过程如下：

1、定义一n维向量U_n；

2、将B’与向量U_n相乘，即：

$B^{\prime }·u_n=A·u_n-\frac{k^T·u_n}{2m}k-\mathrm{diag}(d_1,L,d_N)·u_n$

3、进行迭代计算，即：

$u_{n+1}=\frac{B^{\prime }·u_n}{\left|\right|B^{\prime }·u_n\left|\right|}$

当n→∞时，X≈U_n；β≈||U_n||

式中，X为待求矩阵的绝对值最大特征值对应的特征向量，β为待求矩阵的主特征值，即绝对值最大特征值。

本发明利用图形处理器的并行计算的特性，将网络划分中的特征值和特征向量的计算分解为以下几种基本运算：稀疏矩阵与向量乘法，向量间乘法，向量间加法与数乘，向量间直乘，以及向量范数的计算。通过将计算特征值分解为这几种基本运算，不但将原本庞大的稠密矩阵计算问题简化为稀疏矩阵和向量之间的运算，而且通过基本运算间的合理组合节省了内存空间，并且减少了数据传输的消耗。

以上，仅为本发明的较佳实施例，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求所界定的保护范围为准。