CN103853835A - 基于gpu加速的网络社区检测方法 - Google Patents

Abstract

本发明公开了一种基于GPU加速的网络社区检测方法，主要解决传统NMF网络社区检测方法的运行时间过长，占用空间过多的问题。其实现过程为：(1)构造网络的邻接矩阵；(2)检查网络类型是否满足要求并对邻接矩阵进行预处理；(3)对预处理之后的邻接矩阵进行稀疏表示；(4)对GPU设备进行初始化；(5)将稀疏表示的邻接矩阵传入到GPU设备中；(6)在GPU中进行网络社区检测；(7)将得到的检测结果从GPU设备传到内存中进行归一化处理，得到网络的重叠划分；(8)由重叠划分得到网络的硬化分。本发明利用GPU对NMF社区检测方法进行并行加速，大幅降低了社区检测的运行时间及存储空间，并可处理更大规模的网络数据。

Description

基于GPU加速的网络社区检测方法

技术领域

本发明属于网络技术领域，特别涉及一种网络社区的检测方法，用于复杂网络数据的结构分析。

背景技术

当今社会中，人们的生活总是被各种各样的网络所包围，如社交网、生物网、Web网络等。复杂网络可以用图的形式表示，其中图中的节点表示网络中的数据对象，而图中的边则表示网络对象间的连接情况。一个网络可以用G＝(V,E,A)来表示，其中V是一个包含着n个节点的集合，E是一个包含m条边的集合，A是一个n×n维的邻接矩阵，A_ij代表节点i和节点j之间的边的权重。网络社区检测就是将节点集合V划分到不同的社区中，以得到网络的划分C＝(c₁,c₂,....c_k)，其中k为社区的个数。在重叠网络社区检测中，一个节点可以属于不同的社区，由规模为n×k的隶属度矩阵M来保存其划分的结果，其中m_ij代表了第i个节点属于第j个社区的概率。

网络社区检测对分析网络的拓扑结构、理解复杂网络的功能、发现复杂网络中的隐藏规律和预测复杂网络行为有着重要的理论意义和广泛的应用前景。在社交网站中，可以通过对网络中的用户好友信息对用户进行聚类，可以向用户推荐好友；在电子商务领域，将消费者看作节点，在购买过相同产品的消费者之间建立连边，这样处于同一个社区的消费者可能就具有相同的购物兴趣，进而可以有针对性的向其推荐商品。

非负矩阵分解NMF算法是解决网络社区检测问题的一个经典方法。该算法的思想是把网络社区检测的问题看作一个维数约减的问题。对网络的邻接矩阵A进行分解，分解得到两个矩阵W和H相乘，得到的分解之后的矩阵W即为节点属于不同社区的隶属度矩阵。NMF不仅能够得到网络中比较好的划分结果，同时可以得到每个节点属于每个社区的概率，从而可以同时得到网络的重叠以及硬划分结果。因此NMF算法成为了解决网络社区检测问题的一个经典算法并被广泛使用。

NMF算法存在一个严重的缺点即时间复杂度比较高。经过分析发现其时间复杂度为O(n²×k)，空间复杂度为O(n²)，其中n代表网络中的节点个数，k代表网络的社区数目。当网络规模扩大时，算法运行时间将非常长。因此非常有必要设计并行NMF算法来缩短算法的运行时间。

图形处理器GPU原本是处理计算机图形的专用设备。近十年来，由于高清晰度复杂图形实时处理的要求，GPU发展成为高并行性、多线程、多核的处理器。相比于CPU计算，GPU在进行计算密集型和高度并行性的计算任务上显示了强大的优势。目前，提高CPU性能的主要手段是集成更大的缓存和在同一芯片上集成更多的核，而GPU却可以有效利用晶体管资源，目前主流GPU的运算能力已超过主流通用CPU。传统的基于GPU通用计算的开发方式是将数据打包成纹理，将并行计算任务映射成图形学中的图像着色，使用图形学API进行开发，并利用GPU中的可编程顶点着色单元和像素着色单元作为处理器完成通用计算任务，这个处理过程与图形硬件紧密相关，程序实现非常艰涩并需要大量的底层硬件知识。近年来，利用GPU进行通用计算越来越易于实现，一些标准被相继提出来。

目前的并行程序设计有3种主流框架，分别是CUDA、OpenCL和Direct Compute。本发明采用的是在AMD显卡下的OpenCL平台。OpenCL是由Apple公司提出，并交由Khronos组织进行标准化的开放GPU计算标准语言。OpenCL不仅仅是一种编程语言，更是一个完整的并行编程框架，包括编程语言，API，函数库以及运行时系统来支持软件在整个平台上的开发。在OpenCL编程中，通过核函数来对GPU设备进行编程。核函数运行在GPU上，并完成一个相对独立的功能。由CPU通过运行时函数来调用GPU端的核函数来完成计算。

有学者提出利用GPU设备来对NMF算法进行加速。但是所提出的算法只是对经典NMF模型进行并行。由于网络社区检测问题具有其特殊性：网络邻接矩阵普遍规模比较大而且均为稀疏矩阵而且需要自适应找出划分数目k，套用经典NMF算法往往导致划分结果不佳，加速效果有限，并且不能用于对网络社区的检测。

发明内容

本发明的目的在于针对上述现有方法的不足，提出一种基于GPU加速的网络社区检测方法，以节省计算所需的存储空间，提高社区检测速度。

为实现上述目的，本发明的技术方案包括如下步骤：

(1)构造大小为n×n的网络邻接矩阵A，n为网络中节点的数目；

(2)检查网络是否满足无方向无权重的条件特征，若满足该条件，则去除网络中的孤立节点，即网络中与其他任何节点都没有联系的节点，否则该网络无法处理，退出程序；

(3)对邻接矩阵A采用行格式存储法CSR进行稀疏表示，得到列坐标向量Aj和行首位置向量Ap，同时采用列格式存储法CSC进行稀疏表示，得到行坐标向量Bj和列首位置向量Bp；

(4)利用OpenCL编程框架对GPU设备进行初始化：

选择OpenCL平台，得到OpenCL设备，创建命令队列，创建内存对象，创建OpenCL程序对象，创建核函数对象；

(5)利用0-1之间的随机数分别构造大小为n×k的左分解矩阵W、大小为k×n的右分解矩阵H、大小为1×k的中间向量β，其中k为预设的社区数目，将所述W、H、β、Ap、Aj、Bp和Bj传输到GPU全局内存中，其中，k为初始的网络社区数目；

(6)在GPU中进行网络社区结构检测：

6a)根据上面得到的W、H矩阵以及β向量，计算第一临时变量：T1＝W/(1_n×nH^T+Wdiag(β))，其中1_n×n代表大小为n×n的单位矩阵，H^T代表H矩阵的转置，diag(β)代表β向量的对角矩阵，/代表矩阵间的点除操作；

6b)根据步骤(2)得到的A_p和A_j，计算第二临时变量：T2＝A/(WH)，其中A为网络的邻接矩阵；

6c)根据第二临时变量T2，计算第三临时变量：T3＝T2*H^T；

6d)根据第一临时变量T1和第三临时变量T3对左分解矩阵W进行更新，即

$W\←T1\·T3=\frac{W}{1_{n\×n}{H}^T+\mathrm{Wdiag}\left(\mathrm{\β}\right)}\·\left[\left(\frac{A}{\mathrm{WH}}\right)H^T\right],$

其中·代表两个矩阵的点乘操作，←表示用计算结果替换掉原W矩阵；

6e)计算第四临时变量：T4＝(H·H)1_n×1，其中1_n×1代表大小为n×1的单位向量；

6f)根据步骤6d)得到的左邻接矩阵W，计算第五临时变量：T5＝1_1×n(W·W)，其中1_1×n代表大小为1×n的单位向量；

6g)根据第四临时变量T4及第五临时变量T5对中间向量β中的元素进行更新：β＝(2n+a)/(T4+T5+b)其中，分子参数a固定为8，分母参数b固定为2；

6h)根据左分解矩阵W、右分解矩阵H、中间向量β，计算第六临时变量：T6＝H/(W^T1_n×n+diag(β)H)，W^T代表W矩阵的转置；

6i)根据步骤(2)得到的B_p、B_j，计算第七临时变量：T7＝A/(WH)；

6j)根据第七临时变量T7计算第八临时变量：T8＝W^T*T7；

6k)根据第六临时变量T6和第八临时变量T8，对右分解矩阵H进行更新，即：

$H\←T6\·T8=\frac{W}{{W^{T}1}_{n\×n}+\mathrm{diag}\left(\mathrm{\β}\right)H}\·\left[W^T\left(\frac{A}{\mathrm{WH}}\right)\right],$

6l)重复执行步骤6a)到6k)共100次，得到最终的左分解矩阵W'；

(7)将最终的左分解矩阵W’通过GPU传到计算机内存中；

(8)在计算机内存中对所述左矩阵W'中的元素按行进行归一化，得到网络的重叠划分矩阵M：

$M=\left[\begin{array}{cccc}{m}_{11}& m_{12}& ...& m_{1k}\\ m_{21}& m_{22}& ...& m_{21}\\ .& .& & .\\ .& .& m_{\mathrm{ij}}& .\\ .& .& & .\\ m_{n1}& m_{n2}& ...& m_{\mathrm{nk}}\end{array}\right]$

$m_{\mathrm{ij}}=\frac{{w^{\′}}_{\mathrm{ij}}}{{\mathrm{\Σ}}_{j=1}^k{w^{\′}}_{\mathrm{ij}}}$

其中，w'_ij表示最终的左分解矩阵W'第i行第j列的元素，m_ij代表重叠划分矩阵M的第i行第j列的元素；

(9)根据重叠划分矩阵M得到网络的重叠划分：对于网络中的第i个节点，若M中第i行的任一元素m_ij不等于0，则将该节点划分到第j个社区中且隶属概率为m_ij，否则,该节点不属于第j个社区；

(10)根据网络的重叠划分，将节点划分到隶属概率最大的社区，得到网络的硬化分结果。

本发明具有以下优点：

1，本发明使用非负矩阵分解方法NMF对网络结构进行检测，与现有的网络社区检测技术相比可以同时得到网络的硬划分及重叠划分结果，并可以自适应找到网络的社区数目。

2，本发明采用基于GPU对原非负矩阵分解社区检测方法进行加速，大幅度提升了社区的检测速度。

3，本发明采用行格式存储法CSR与列格式存储法CSC对网络的邻接矩阵进行稀疏表示，节省了处理过程所需的存储空间，使得本发明可处理更大规模的网络数据。

附图说明

图1是本发明的流程图；

图2是用本发明处理两个真实网络的结果图；

图3是用本发明处理科学家协作网络的加速效果示意图。

具体实施方式

下面结合附图，对本发明做进一步的描述。

参照图1，本发明具体实施步骤如下：

步骤1：构造网络邻接矩阵A。

邻接矩阵A的大小为n×n的，n为网络中节点的数目；

若网络中的节点i与节点j有连接时，邻接矩阵A中的第i行第j列的元素a_ij=1；

若节点i与节点j无连接，则a_ij=0。

步骤2：检查网络是否满足无方向无权重的条件特征。

遍历邻接矩阵A中的所有元素，若

a_ij＝a_ji则该网络满足无方向特征，若

a_ij∈{0,1}，则该网络满足无权重特征；

若网络同时满足无方向无权重条件，则去除网络中的孤立节点，否则该网络无法处理，退出程序。

网络中的孤立节点是指与其他任何节点都没有联系的节点，若第i个节点是孤立节点，去除该孤立节点是指删除邻接矩阵A中的第i行和第i列的所有元素。

步骤3：对邻接矩阵A采用行格式存储法CSR进行稀疏表示，得到列坐标向量Aj和行首位置向量Ap，同时采用列格式存储法CSC进行稀疏表示，得到行坐标向量Bj和列首位置向量Bp。

3a)将邻接矩阵A的所有非零元素所在的列坐标，按照行优先的顺序保存为列坐标向量Aj；

3b)将邻接矩阵A每行第一个非零元素在Aj中的位置保存到行首位置向量Ap中，并将加入到网络的节点数目作为行首位置向量Ap的最后一个元素；

3c)将邻接矩阵A的所有非零元素所在的行坐标按照列优先的顺序保存为Bj；

3d)将邻接矩阵A每列第一个非零元素在Bj中的位置保存到列首位置向量Bp，并将加入到网络的节点数目作为列首位置向量Bp的最后一个元素。

步骤4：利用OpenCL编程框架对GPU设备进行初始化。

选择OpenCL平台，得到OpenCL设备，创建命令队列，创建内存对象，创建OpenCL程序对象，创建核函数对象。

这一系列操作为OpenCL标准所要求，依次调用OpenCL提供的底层运行时函数来完成。

步骤5：数据初始化。

利用0-1之间的随机数分别构造左分解矩阵W、右分解矩阵H和中间向量β，其中左分解矩阵W的大小为n×k，右分解矩阵H的大小为k×n，中间向量β的大小为1×k，k为预设的社区数目；

将所述左分解矩阵W、右分解矩阵H、中间向量β、列坐标向量Aj、行首位置向量Ap、列首位置向量Bj和列首位置向量Bp传输到GPU的全局内存中。

步骤6：在GPU中进行网络社区检测。

网络社区检测的过程即为对左分解矩阵W、右分解矩阵H、中间向量β多次更新迭代更新的过程，一次迭代需要执行的计算为：

$W\←T1\·T3=\frac{W}{1_{n\×n}{H}^T+\mathrm{Wdiag}\left(\mathrm{\β}\right)}\·\left[\left(\frac{A}{\mathrm{WH}}\right)H^T\right]$

$H\←T6\·T8=\frac{W}{{W^{T}1}_{n\×n}+\mathrm{diag}\left(\mathrm{\β}\right)H}\·\left[W^T\left(\frac{A}{\mathrm{WH}}\right)\right]$

$\mathrm{\β}\←\frac{2n+a}{1_{n\×n}(W\·W)+(H\·H)1_{n\×1}+b},$

其中，1_n×n代表大小为n×n的单位矩阵，H^T代表分解右矩阵H的转置，diag(β)代表中间向量β的对角矩阵，·代表两个矩阵的点乘操作，←表示用计算结果对左边变量进行更新，W^T代表分解右矩阵W的转置，1_1×n代表大小为1×n的单位向量，1_n×1代表大小为n×1的单位向量，分子参数a固定为8，分母参数b固定为2。

对上述三个公式进行拆解，在GPU中完成各拆解之后的运算，其执行步骤如下：

6a)计算第一临时变量：T1＝W/(1_n×nH^T+Wdiag(β))，本步骤在GPU中的实现方式如下：

6a1)在GPU中开辟n×k个线程，其中n为网络节点的数目，k为初始的社区数目，每个线程负责计算第一临时变量T1中的一个元素；

6a2)在上面开辟的线程中，对于位置为第x行第y列的线程，计算第一临时变量T1中第x行第y列的元素:并将计算结果T1_xy写到第一临时变量T1的第x行第y列的位置上，其中h_yi代表右分解矩阵H矩阵第y行第i列的元素，β_x代表中间向量β的第x个元素，w_xy代表左分解矩阵W的第x行第y列的元素；

6b)计算第二临时变量：T2＝A/(WH)，本步骤在GPU中的实现方式如下：

6b1)在GPU中开辟n×n个线程；

6b2)在步骤6b1)开辟的线程中，对于位置为第x行第y列的线程，计算其负责处理的元素在列坐标向量A_j中的位置：index＝Ap_x+y，及该元素在第二临时变量T2中的行坐标row=x、列坐标col＝Aj_index，其中Ap_x代表行首位置向量Ap的第x个元素，Aj_index代表列坐标向量Aj的第index个元素；

6b3)读取左分解矩阵W第row行的元素，并放到局部内存As中；读取右分解矩阵H中第col列的元素，并放到局部内存Bs中；对这两个内存As与Bs中保存的向量求内积，并对所求内积结果取倒数，将取倒数后的最终结果写到第二临时变量T2的第index个位置上；

6c)计算第三临时变量：T3＝T2*H^T，本步骤在GPU中的实现方式如下：

6c1)在GPU中开辟n×k个线程；

6c2)在步骤6c1)开辟的线程中，对于位于第x行第y列的线程，计算邻接矩阵A的第x行中非零元素的个数zn＝Ap_x+1-Ap_x，读取第二临时变量T2中相应的元素并放入局部内存As中，即

i＝1,2,...,zn，其中As[i]表示局部内存As第i个位置保存的元素值，表示第二临时变量T2中第i+Ap_x个元素；将分解右矩阵H的第y行元素放入局部内存Bs中，计算第三临时变量T3第x行第y列的元素：

${T3}_{\mathrm{xy}}=\underset{i=1}{\overset{\mathrm{zn}}{\mathrm{\Σ}}}\mathrm{As}\left[i\right]*\mathrm{Bs}\left[{\mathrm{Aj}}_x\right]$

其中，Bs[Aj_x]代表局部内存Bs第Aj_x个位置保存的元素值，Aj_x表示列坐标向量Aj的第x个元素值；

6d)用第一临时变量T1和第三临时变量T3更新左分解矩阵W：

$W\←T1\·T3=\frac{W}{1_{n\×n}{H}^T+\mathrm{Wdiag}\left(\mathrm{\β}\right)}\·\left[\left(\frac{A}{\mathrm{WH}}\right)H^T\right],$

本步骤在GPU中的实现如下：

6d1)在GPU中开辟n×k个线程；

6d2)在步骤6d1)开辟的线程中，对于位置为第x行第y列的线程，计算左分解矩阵W中第x行第y列元素的更新值：

其中T1_xy代表第一临时变量T1的第x行第y列的元素，T3_xy代表第三临时变量T3的第x行第y列的元素，将计算结果写到左分解矩阵W第x行第y列的位置上；

6e)计算第四临时变量：T4＝(H·H)1_n×1，本步骤在GPU中的实现如下：

6e1)在GPU中开辟k×1个线程；

6e2)在步骤6e1)开辟的线程中，对于第x行的线程，计算第四临时变量T4的第x个元素：

并将计算结果写到第四临时变量T4第x个位置上，其中h_xy代表右分解矩阵H第x行第y列的元素；

6f)计算第五临时变量：T5＝1_1×n(W·W)，本步骤在GPU中的实现如下：

6f1)在GPU中开辟1×k个线程；

6f2)在步骤6f1)开辟的线程中，对于第y列的线程，计算第五临时变量T5的第y个元素：

并将结果写到第五临时变量T5的第y个位置上，其中w_xy代表左分解矩阵W第x行第y列的元素；

6g)根据第四临时变量T4及第五临时变量T5对中间向量β进行更新，本步骤在GPU中的实现如下：

6g1)在GPU中开辟k×1个线程；

6g2)在步骤6g1)开辟的线程中，对于第x列的线程，计算中间向量β的第x个元素的更新值：其中T4_x代表第四临时变量T4的第x个元素，T5_x代表第五临时变量T5的第x个元素，将结果写到中间变量β的第x个位置上；

6h)计算第六临时变量：T6＝H/(W^T1_n×n+diag(β)H)，本步骤在GPU中的实现如下：

6h1)在GPU中开辟k×n个线程；

6h2)在6h1)中开辟的线程中，对于位置为第x行第y列的线程，计算第六临时变量T6中第x行第y列的元素：

并将计算结果T6_xy写到第六临时变量T6的第x行第y列的位置上，其中β_y代表中间变量β的第y个元素；

6i)计算第七临时变量：T7＝A/(WH)，本步骤在GPU中的实现如下：

6i1)在GPU中开辟n×n个线程；

6i2)在步骤6i1)开辟的线程中，对于位置为第x行第y列的线程，计算其负责处理的元素在行坐标向量Bj中的位置：Index＝Bp_y+x，及该元素在第七临时变量T7中的行坐标Row＝Bj_Index、列坐标Col＝y，其中Bp_y代表列首位置向量Bp的第y个元素，Bj_Index代表行坐标向量Bj的第Index个元素；

6i3)读取左分解矩阵W第Row行的元素，并放到局部内存As中；读取右分解矩阵H中第Col列的元素，并放到局部内存Bs中；对这两个内存As与Bs中保存的向量求内积，并对所求内积结果取倒数，将取倒数后的最终结果写到第七临时变量T7的第Index个位置上；

6j)计算第八临时变量：T8＝W^T*T7，本步骤在GPU中的实现如下：

6j1)在GPU中开辟k×n个线程；

6j2)在步骤6j1)开辟的线程中，对于位于第x行第y列的线程，计算邻接矩阵A的第y列中非零元素的个数Zn＝Bp_y+1-Bp_y，读取左分解矩阵W第x列的元素放入局部内存As中，读取第七临时变量T7中的相应元素并放入局部内存Bs中，即

i＝1,2,...,Zn,其中Bs[i]表示局部内存Bs第i个位置保存的元素值，

表示第七临时变量T7中的第i+Bp_y个元素；计算第八临时变量T8第x行第y列的元素：

${T8}_{\mathrm{xy}}=\underset{i=1}{\overset{\mathrm{zn}}{\mathrm{\Σ}}}\mathrm{As}\left[{\mathrm{Bj}}_y\right]*\mathrm{Bs}\left[i\right]$

其中As[Bj_y]代表局部内存As中第Bj_y个位置保存的元素值，Bj_y表示行坐标向量Bj的第y个元素；

6k)用第六临时变量T6和第八临时变量T8更新右分解矩阵H：

$H\←T6\·T8=\frac{W}{{W^{T}1}_{n\×n}+\mathrm{diag}\left(\mathrm{\β}\right)H}\·\left[W^T\left(\frac{A}{\mathrm{WH}}\right)\right],$

本步骤在GPU中的实现如下：

6k1)在GPU中开辟n×k个线程；

6k2)在步骤6k1)开辟的线程中，对于位置为第x行第y列的线程，计算左分解矩阵H中第x行第y列元素的更新值：

其中T6_xy代表第六临时变量T6的第x行第y列的元素，T8_xy代表第三临时变量T8的第x行第y列的元素，将计算结果写到左分解矩阵H第x行第y列的位置上；

6l)重复执行步骤6a)到6k)共100次，得到最终的左分解矩阵W'；

步骤7：将最终的左分解矩阵W’通过GPU传到计算机内存中,在计算机内存中对该最终的左分解矩阵W'中的元素按行进行归一化，得到网络的重叠划分矩阵M：

$M=\left[\begin{array}{cccc}{m}_{11}& m_{12}& ...& m_{1k}\\ m_{21}& m_{22}& ...& m_{21}\\ .& .& & .\\ .& .& m_{\mathrm{ij}}& .\\ .& .& & .\\ m_{n1}& m_{n2}& ...& m_{\mathrm{nk}}\end{array}\right]$

$m_{\mathrm{ij}}=\frac{{w^{\′}}_{\mathrm{ij}}}{{\mathrm{\Σ}}_{j=1}^k{w^{\′}}_{\mathrm{ij}}}$

其中，w'_ij表示最终的左分解矩阵W'第i行第j列的元素，m_ij代表重叠划分矩阵M的第i行第j列的元素。

步骤8：根据重叠划分矩阵M得到网络的重叠划分。

对于网络中的第i个节点，若M中第i行的任一元素m_ij不等于0，则将该节点划分到第j个社区中且隶属概率为m_ij，否则，该节点不属于第j个社区。

步骤9：根据网络的重叠划分，将节点划分到隶属概率最大的社区，得到网络的硬化分结果。

本发明的效果可通过以下仿真进一步说明：

1仿真内容：应用本发明方法对12种不同的网络进行处理，并使用原始NMF算法进行对比。测试本发明方法得到的网络社区结构的有效性，以及其相对于原始NMF算法的加速效果。

2.仿真条件：

本发明的仿真硬件环境：CPU为AMD4Dual Core2.8GHz，内存为32GB DDR2，GPU为AMD Radeon HD77701GHz显存为1GB DDR5；

本发明的仿真软件环境：编译环境为Microsoft Visual Studio2010,GPU端代码标准为OpenCL1.2版本，操作系统为Windows764位系统。

3仿真实验结果

仿真1.本发明方法进行网络社区检测，测试其有效性。

将本发明方法应用于12种不同规模的网络中，用评价指标NMI和模块度Q来测试所得结果的准确性。结果在表1中给出，左边三列中是人工的网络数据的结果，右列两列为真实网络数据的结果。从表中的结果可以看出算法对于人工网络和真实网络均可以得到比较好的结果。图2给出其中两个网络的具体划分结果，其中图2(a)是Netscience网络的一部分划分结果，图2(b)是LFR500网络的社区结构，在图2中用虚线椭圆框来标记划分结果，两个节点在不同框内代表它们被划分到不同的社区。从该图可以看出，本发明方法可以得到比较合理的网络社区划分。

表1.采用发明方法对处理12种网络的结果

人工网络	Q	NMI	真实网络	Q
					Benchmark	0.4501	1	Facebook	0.62
LFR(500)	0.8131	1	Email	0.5312
					LFR(1000)	0.9039	0.9714	NetScience	0.905
LFR(5000)	0.9075	0.9435	Power	0.599
					LFR(10000)	0.8599	0.9329	Scientists	0.647
LFR(50000)	0.8999	0.9191	Hep	0.7716

仿真2.测试本发明方法相比于原始NMF方法的加速效果

在表2中统计了原始NMF算法和本发明方法在处理上述12种数据花费的时间，其中n代表网络的节点数目，k代表初始的社区数目，T_NMF代表原NMF方法处理该网络需要的时间，T_SNMF代表本发明方法处理该网络需要的时间，时间单位均为秒，加速比speedup ratio等于T_NMF/T_SNMF表示本发明方法相对于原NMF方法的加速性能。表2中的*号表示由于运行时间以及存储空间的限制，该算法对该网络无法处理，也无法计算加速比。

表2.采用发明方法对处理8种网络的加速效果

由表2可知，本发明方法取得了较高的加速比，如对于节点数目为5000的LFR网络，原始NMF方法的运行时间超过8个小时，本发明方法运行时间仅为71秒。原始NMF方法对节点数5000以上的网络便很难进行处理，而本发明可以处理超过50000万个节点的网络。另外，从表格的最后一列可以看出本发明方法的加速比会随着网络规模的增大而提高。

仿真3.测试初始社区数目k对本发明方法加速性能的影响

本仿真选择Netscience网络为例进行测试，图3给出了本发明方法在不同的k值下所达到的加速比。从该图看出本发明方法的加速比随着k值的增大而增大。k取256时，加速比已经达到了381，远大于k取64时的201。因此，利用本发明方法可以有效缓解k值增大对算法运行时间的影响。在实际网络划分中，社区数目普遍比较高，这一效果会更加明显。

Claims (4)

1.一种基于GPU加速的网络社区检测方法，包括如下步骤：

(1)构造大小为n×n的网络邻接矩阵A，n为网络中节点的数目；

(2)检查网络是否满足无方向无权重的条件特征，若满足该条件，则去除网络中的孤立节点，即网络中与其他任何节点都没有联系的节点，否则该网络无法处理，退出程序；

(3)对邻接矩阵A采用行格式存储法CSR进行稀疏表示，得到列坐标向量Aj和行首位置向量Ap，同时采用列格式存储法CSC进行稀疏表示，得到行坐标向量Bj和列首位置向量Bp；

(4)利用OpenCL编程框架对GPU设备进行初始化：

选择OpenCL平台，得到OpenCL设备，创建命令队列，创建内存对象，创建OpenCL程序对象，创建核函数对象；

(5)利用0-1之间的随机数分别构造大小为n×k的左分解矩阵W、大小为k×n的右分解矩阵H、大小为1×k的中间向量β，其中k为预设的社区数目，将所述W、H、β、Ap、Aj、Bp和Bj传输到GPU全局内存中；

(6)在GPU中进行网络社区结构检测：

6a)根据上面得到的W、H矩阵以及β向量，计算第一临时变量：T1＝W/(1_n×nH^T+Wdiag(β))，其中1_n×n代表大小为n×n的单位矩阵，H^T代表H矩阵的转置，diag(β)代表β向量的对角矩阵，/代表矩阵间的点除操作；

6b)根据步骤(2)得到的A_p和A_j，计算第二临时变量：T2＝A/(WH)，其中A为网络的邻接矩阵；

6c)根据第二临时变量T2，计算第三临时变量：T3＝T2*H^T；

6d)根据第一临时变量T1和第三临时变量T3对左分解矩阵W进行更新，

即 $W\←T1\·T3=\frac{W}{1_{n\×n}{H}^T+\mathrm{Wdiag}\left(\mathrm{\β}\right)}\·\left[\left(\frac{A}{\mathrm{WH}}\right)H^T\right],$

其中·代表两个矩阵的点乘操作，←表示用计算结果替换掉原W矩阵；

6e)计算第四临时变量：T4＝(H·H)1_n×1，其中1_n×1代表大小为n×1的单位向量；

6f)根据步骤6d)得到的左邻接矩阵W，计算第五临时变量：T5＝1_1×n(W·W)，其中1_1×n代表大小为1×n的单位向量；

6g)根据第四临时变量T4及第五临时变量T5对中间向量β进行更新：β＝(2n+a)/(T4+T5+b)

其中，分子参数a固定为8，分母参数b固定为2；

6h)根据左分解矩阵W、右分解矩阵H、中间向量β，计算第六临时变量：T6＝H/(W^T1_n×n+diag(β)H)，W^T代表W矩阵的转置；

6i)根据步骤(2)得到的B_p、B_j，计算第七临时变量：T7＝A/(WH)；

6j)根据第七临时变量T7计算第八临时变量：T8＝W^T*T7；

6k)根据第六临时变量T6和第八临时变量T8，对右分解矩阵H进行更新，即：

$H\←T6\·T8=\frac{W}{{W^{T}1}_{n\×n}+\mathrm{diag}\left(\mathrm{\β}\right)H}\·\left[W^T\left(\frac{A}{\mathrm{WH}}\right)\right],$

6l)重复执行步骤6a)到6k)共100次，得到最终的左分解矩阵W'；

(7)将最终的左分解矩阵W’通过GPU传到计算机内存中，对该最终的左分解矩阵W'中元素按行进行归一化，得到网络的重叠划分矩阵M：

$M=\left[\begin{array}{cccc}{m}_{11}& m_{12}& ...& m_{1k}\\ m_{21}& m_{22}& ...& m_{21}\\ .& .& & .\\ .& .& m_{\mathrm{ij}}& .\\ .& .& & .\\ m_{n1}& m_{n2}& ...& m_{\mathrm{nk}}\end{array}\right]$

$m_{\mathrm{ij}}=\frac{{w^{\′}}_{\mathrm{ij}}}{{\mathrm{\Σ}}_{j=1}^k{w^{\′}}_{\mathrm{ij}}}$

其中，w'_ij表示最终的左分解矩阵W'第i行第j列的元素，m_ij代表重叠划分矩阵M的第i行第j列的元素；

(8)根据重叠划分矩阵M得到网络的重叠划分：对于网络中的第i个节点，若M中第i行的任一元素m_ij不等于0，则将该节点划分到第j个社区中且隶属概率为m_ij，否则,该节点不属于第j个社区；

(9)根据网络的重叠划分，将节点划分到隶属概率最大的社区，得到网络的硬化分结果。

2.根据权利要求书1所述的基于GPU加速的网络社区检测方法，其特征在于，步骤6a)所述的计算第一临时变量T1操作，按如下步骤进行：

6a1)在GPU中开辟n×k个线程，其中n为网络节点的数目，k为初始的社区数目，每个线程负责计算T1中的一个元素；

6a2)在上面开辟的线程中，任意选取其中第x行第y列的线程来计算T1中第x行第y列的元素:

并将计算结果T1_xy写到T1变量的第x行第y列的位置上，其中h_yi代表右分解矩阵H矩阵第y行第i列的元素，β_x代表中间向量β的第x个元素，w_xy代表左分解矩阵W的第x行第y列的元素。

3.根据权利要求书1所述的基于GPU加速的网络社区检测方法，其特征在于，步骤6b)所述的计算T2操作，按如下步骤进行：

6b1)在GPU中开辟n×n个线程；

6b2)在上面开辟的线程中，对于位置为第x行第y列的线程，计算其负责处理的元素在列坐标向量A_j中的位置：index＝Ap_x+y,及该元素在结果矩阵中的行坐标row=x、列坐标col＝Aj_index，其中Ap_x代表行首位置向量Ap的第x个元素，Aj_index代表列坐标向量Aj的第index个元素；

6b3)读取左分解矩阵W第row行的元素，并放到局部内存As中，读取右分解矩阵H中第col列的元素，并放到局部内存Bs中，对As与Bs中保存的向量求内积，对所求内积结果取倒数，将最终结果写到第二临时变量T2的第index个位置上。

4.根据权利要求书1所述的基于GPU加速的网络社区检测方法，其特征在于，步骤6c)所述的计算第三临时变量T3，按如下步骤进行：

6c1)在GPU中开辟n×k个线程；

6c2)在上面开辟的线程中，对于位于第x行第y列的线程，计算邻接矩阵A的第x行中非零元素的个数zn＝Ap_x+1-Ap_x，读取第二临时变量T2中相应的元素并放入局部内存As中，即

i＝1,2,...,zn，其中As[i]表示局部内存As第i个位置保存的元素值，

表示第二临时变量T2中第i+Ap_x个元素，将H矩阵的第y行元素放入局部内存Bs中，计算第三临时变量T3第x行第y列的元素：

${T3}_{\mathrm{xy}}=\underset{i=1}{\overset{\mathrm{zn}}{\mathrm{\Σ}}}\mathrm{As}\left[i\right]*\mathrm{Bs}\left[{\mathrm{Aj}}_x\right]$

其中，Bs[Aj_x]代表Bs第Aj_x个位置保存的元素值，Aj_x表示列坐标向量Aj的第x个元素值。

Images