LabelRank非重叠社区发现算法介绍及代码实现（A Stabilized Label Propagation Algorithm for Community Detection in Networks）

最近在研究基于标签传播的社区分类，LabelRank算法基于标签传播和马尔科夫随机游走思路上改装的算法，引用率较高，打算将代码实现，便于加深理解。

这个算法和Label Propagation 算法不同的是计算复杂度较高，对每个标签都确定了概率，但是准确性比Label Propagation算法好。

一、概念

相关概念不再累述，详情见前两篇文章

二、算法思路

首先建立一个标签集合，C={1，2，……n}，n是节点的数量。标签概率向量P_i（1*n），P_i(c)=节点i对标签c的概率估计，迭代过程中每个节点的对标签c的概率估计等于其邻居节点对标签c的概率估计平均，详见公式（1）

有此可得n*n维标签概率矩阵P（i→j）=[p_1,p_2,...p_n],迭代过程可以用矩阵乘法表示A*P，其中A是网络的邻接矩阵（01矩阵）。这个思路其实可以追溯到eigenvector Centrality算法¹，文献1已证明P会收敛下来。就这样就完了吗？并没有看到如何传递标签或者选择标签？

作者做的就是不停地缩放P中元素，然后删除一些概率较小的标签从P中，不停地减少标签个数，知道每个节点的标签序列不再变化，迭代停止，拥有最大概率的标签就是节点所属的社区。具体流程见下

（1）Propagation

初始阶段，每个节点访问邻居概率皆相等，见公式（3），每次迭代即左乘上一阶段的P，得到本阶段节点对每个标签的预估概率

（2）Inflation

根据公式（2）不停地迭代，矩阵中0，计算复元素逐渐被取代，复杂度越来越高，流程（2）和（3）就是为减少复杂度而做的工作。首先利用公式（4）将矩阵中的元素极端处理，使值大的越来越大，值小的越来越小

（3）Cut off

这一阶段就是在公式（4）的基础上进行删除操作，将P中低于r的阈值全都置换成0，最终得到的P参与下一次迭代

（4）Explicit Conditional Update

减少算法的另一个途径就是满足某一条件的节点停止更新，具体操作就是如果节点的最大标签（对n个标签估计概率最高的那个标签）和他的邻居节点最大标签的吻合度高于q（提前给出，一般去0.7左右），那么这个节点就可以停止更新了

（5）Stop Criterion

每个节点的最大评估概率的标签不再变化，迭代停止，具有相同标签的节点归为一个社区

三、参考文献

[1]Poulin R, Boily M C, Mâsse B R. Dynamical systems to define centrality in social networks[J]. Social Networks, 2000, 22(3):187-220.

Dynamical systems to define centrality in social networks

[2]Xie J, Szymanski B K. LabelRank: A stabilized label propagation algorithm for community detection in networks[C]// Network Science Workshop. IEEE, 2013:138-143.

A Stabilized Label Propagation Algorithm for Community Detection in Networks

四、代码（matlab）

代码目前还有一点点问题，后期调试后再更新

function [R,count]=LabelR(A,in,r,q)
%  LabelRank LabelRank: " A Stabilized Label Propagation
%  Algorithm for Community Detection in Networks "
%  Author: YY
%  Created on 2017.05.09
% Inputs : 
% A : adjacent matrix 
% in : Inflation parameter
%   : default =2
% q : Conditional Update parameter
%        default = 0.7
% r   : Cut off parameter
%  : default = 0.1
% Output :
% R : community classfication
%% 
   %  Step1 : Propagation
Aori=A;             
A=A+eye(length(A));%  add selfloop
k=repmat(sum(A,2),[1,length(A)]);
P0=A./k;
Ppre=A*P0;
a=1;
COM={};
count=0;
%% 
   % Step2: Inflation
   while a
       Pnow=A*Ppre;
       Pin=Pnow.^in ;
       k=repmat(sum(Pin,2),[1,length(A)]);
       Pnow=Pin./k;
 %% 
   % Step3: Cutoff
       index= Pnow<r;
       Pnow(index)=0;
%% 
   % Step4: Explicit Conditional Update
    MaNow=max(Pnow,[],2);
    MaPre=max(Ppre,[],2);
    restart=[];
       for i=1:length(A)
           gain=0;
           Nb=find( Aori(i,:));
           MaxI=max(Pnow(i,:));
           MaxI=find(Pnow(i,:)==MaxI);
           MaxNb=MaNow(Nb);
           for k=1:length(Nb)
               MaxNbID=find(Pnow(Nb(k),:)==MaxNb(k));
               if all(ismember(MaxI,MaxNbID));% 1,2和1；1和1,2;1,2和1,2,4或者1,3,4
                   gain=gain+1;
               end
           end
           if gain>=q*length(Nb)
              restart=[i,restart]; 
           end
       end
        Pnow(restart,:)=Ppre(restart,:);
 %% 
   % Step5: Stop Criterion
       if all(ismember(find(Pnow(i,:)==MaNow(i)),find(Ppre(i,:)==MaPre(i))))
           a=0;
       end
       Ppre=Pnow;
       count=count+1;
   end
   R=Pnow;
end