Matlab实现小世界网络生成及其分析

引言

本文的主要工作包括：

• 对社会网络的结构属性做了一个介绍
• 对群聚系数、平均路径长度以及节点度分布计算方式做了一个介绍
• 提供了Watts­-Strogatz“小世界”网络生成算法流程及其matlab代码实现

社会网络的结构属性

我们可以用图论的语言来刻画社会网络。此时，一个社会网络可以表示为由点集$V(G)$和边集$E(G)$所共同决定的一个图$G(V,E)$。其中，集合$V(G)$的元素个数表示网络节点的个数或社会群体的规模。
进一步，为了描述不同网络结构的图的表示的特点，学者们提出并研究了许多图的属性特征。常见的图结构属性包括：

1. 节点度：指与该节点相连的边的数目。在无向图中，节点度即为与该节点相连接的所有其他节点的数目。而在有向图中，将所有指向某节点的边的数量记为该节点的入度，所有从该节点出发指向别的节点的边的数量记为该节点的出度。网络平均度反应了网络的疏密程度。
2. 度分布：指随机选定一个节点，该节点的度数等于$k$的概率$p(k)$，即网络中节点的度数等于$k$的比例。通过度分布则可以刻画不同节点的重要性。
3. 群聚系数：用于描述网络中与同一节点相连的节点间也互为相邻节点的程度。其用于刻画社会网络中一个人朋友们之间也互相是朋友的概率，反映了社会网络中的聚集性。
4. 介数：为图中某节点承载整个图所有最短路径的数量，通常用来评价节点的重要程度，比如在连接不同社群之间的中介节点的介数相对于其他节点来说会非常大，也体现了其在社会网络信息传递中的重要程度。
5. 网络直径：指网络中任意两个节点之间的最短路径的最大值，影响着复杂系统的收敛时间和整个网络的稳定性。
6. 平均路径长度：指网络中任意两个节点之间的最短路径的平均值。用来度量整体网络节点之间通信的有效性，即一个个体的信息平均传播多少距离才能被另一个个体接收。

其中，度分布、群聚系数以及平均路径长度是衡量一个网络拓扑结构鲁棒性的三个最重要的指标。

Watts­-Strogatz小世界网络生成算法

前面，我们已经了解到不同的社会网络会具有不同网络结构属性。并且，我们可以根据群聚系数、平均路径长度的差异将网络结构划分为三类：

1. 规则网络：具有明显的局部群聚结构和较大的网络群聚系数，但网络平均路径长度会随网络节点数线性增长。
2. 随机网络：网络平均路径长度随网络节点数对数增长，但不具有局部群聚结构，即网络群聚系数很小。
3. “小世界”网络：具有同规则网络相当的群聚系数以及和随机网络一样较小的网络平均路径长度。

顾名思义，“小世界”网络具有着“小世界现象”，亦即“六度空间理论”。现实社会中的许多网络结构都具有这一特性。例如：蠕虫线虫的神经网络、电影演员的合作关系网络、美国西部的电力网络[1]。

基于此，Watts和Strogatz提出了一种通过重连规则网络中的已有边从而生成小世界网络的算法。
具体的算法流程图如下：

matlab代码实现

%% Copyright 2015 The MathWorks, Inc.
% H = WattsStrogatz(N,K,beta) returns a Watts-Strogatz model graph with N
% nodes, N*K edges, mean node degree 2*K, and rewiring probability beta.
%
% beta = 0 is a ring lattice, and beta = 1 is a random graph.
% Inputs:
% N: total network nodes
% K: node degree
% beta: rewiring probability beta
% Outputs:
% retG: A graph object
% retMat: The adjacency matrix
% Examples:
% [retG,retMat] = WattsStrogatz(N,K,beta); % accept two results
% retG = WattsStrogatz(N,K,beta); % only accept the graph object
% [~, retMat] = WattsStrogatz(N,K,beta); % only accept the adjacency matrix
function [retG, retMat]= WattsStrogatz(N,K,beta)
hwait = waitbar(0,'Please wait. Generating WattsStrogatz graph.');
% Connect each node to its K next and previous neighbors. This constructs
% indices for a ring lattice.
s = repelem((1:N)',1,K);
t = s + repmat(1:K,N,1);
t = mod(t-1,N)+1;

% Rewire the target node of each edge with probability beta
for source=1:N    
    waitbar(source/N,hwait,sprintf('Please wait. Generating WattsStrogatz graph.\n%.1f %%',(source/N)*100));
    switchEdge = rand(K, 1) < beta;
    
    newTargets = rand(N, 1);
    newTargets(source) = 0;
    newTargets(s(t==source)) = 0;
    newTargets(t(source, ~switchEdge)) = 0;
    
    [~, ind] = sort(newTargets, 'descend');
    t(source, switchEdge) = ind(1:nnz(switchEdge));
end

retG = graph(s,t); % 返回图
retMat = adjacency(retG); % 返回邻接矩阵
delete(hwait);
end

网络结构分析

在这一重连过程中，网络结构会随着边重连概率$p$的值的不同而有所改变。

相应地，网络结构属性也会有所改变。我们只考虑网络结构的群聚系数以及平均路径长度的改变情况。并且，统计了不同的边重连概率$p$下的结构属性的变化情况。

图2：群聚系数与平均路径长度随边重连概率变化的曲线图

从图2中，我们可以发现：

1. 当边重连概率$p$在区间[0,0.01]增长时，网络结构的群聚系数不会发生太大的变化，而网络的平均路径长度却会快速下降。这就是“捷径边”的影响。影响了网络的全局结构属性
2. 而当边重连概率进一步增加时，原来网络结构中的边随机性增强。最终在$p=1$的时候得到一个完全随机的网络。在这一过程中，群体局部的群聚结构被打破，进而导致群聚系数逐渐减小。

网络结构属性分析及实现代码

关于群聚系数和平均路径长度的定义式，已经有许多公开的不错文献。大家可以自行参考。以下只贴代码：
计算群聚系数的代码

function [acc, c] = avgClusteringCoefficient(graph)
%Computes the Average Clustering Coefficient for the undirected, unweighted
%graph input as adjacency matrix 'graph' and returns it in variable 'acc'. 
%It also returns the local clustering coefficient of each node in 'c'.
%This implementation does not take into account nodes with zero degree.
%
%The definition of clustering coefficient used in this implementation was
%taken from:
%
%   Watts,D.J. and Strogatz,S.H. (1998) Collective dynamics of 
%   "small-world" networks. Nature, 393, 440-442.
%
%INPUT
%   graph -> The adjacency matrix representation of a graph. It has to be a
%            NxN matrix where N is the number of nodes in the graph. This
%            parameter is required.
%
%OUTPUT
%   acc -> Average clustering coefficient of the input graph.
%   c -> Local clustering coefficient of each of the graph's nodes
%
%Example usage:
%
%   [acc, c] = avgClusteringCoefficient(my_net);
%   acc = avgClusteringCoefficient(my_net);
%   [~, c] = avgClusteringCoefficient(my_net);
%
% Copyright (C) Gregorio Alanis-Lobato, 2014

%% Input parsing and validation
ip = inputParser;

%Function handle to make sure the matrix is symmetric
issymmetric = @(x) all(all(x == x.'));

addRequired(ip, 'graph', @(x) isnumeric(x) && issymmetric(x));

parse(ip, graph);

%Validated parameter values
graph = ip.Results.graph;


%% Clustering coefficient computation

%Make sure the graph unweighted!!!
graph(graph ~= 0) = 1; 

deg = sum(graph, 2); %Determine node degrees
cn = diag(graph*triu(graph)*graph); %Number of triangles for each node

%The local clustering coefficient of each node
c = zeros(size(deg));
c(deg > 1) = 2 * cn(deg > 1) ./ (deg(deg > 1).*(deg(deg > 1) - 1)); 

%Average clustering coefficient of the graph
acc = mean(c(deg > 1)); 
end

计算平均路径长度的代码

function [APL, L] = avgPathLength(adjMat)
%Computes the Average path length for the undirected, unweighted
%graph input as adjacency matrix 'adjMat' and returns it in variable 'APL'. 
%It also returns the global path length of each node in 'L'.
%This implementation does not take into account nodes with zero degree.
%
% Average path length is a concept in network topology that is defined as 
% the average number of steps along the shortest paths for all possible 
% pairs of network nodes. It is a measure of the efficiency of information 
% or mass transport on a network.
% L_g = sum(d(vi,vj))/n*(n-1)
% where d(vi,vj) denotes the shortest length between two nodes, and equals
% to 0 if node i cannot be reached from j

%INPUT
%   adjMat -> The adjacency matrix representation of a graph. It has to be a
%            NxN matrix where N is the number of nodes in the graph. This
%            parameter is required.
%
%OUTPUT
%   APL -> Average Path Length of the input graph.
%   L -> Global Path length of each of the graph's nodes
%
%Example usage:
%
%   [APL, L] = avgClusteringCoefficient(my_net);
%   APL = avgClusteringCoefficient(my_net);
%   [~, L] = avgClusteringCoefficient(my_net);
%
% Copyright (C) KunTang, 2022

%% Input parsing and validation
ip = inputParser;

%Function handle to make sure the matrix is symmetric
issymmetric = @(x) all(all(x == x.'));

addRequired(ip, 'graph', @(x) isnumeric(x) && issymmetric(x));

parse(ip, adjMat);

%Validated parameter values
adjMat = ip.Results.graph;
%% Path Length computation

%Make sure the graph unweighted!!!
adjMat(adjMat ~= 0) = 1; 

graphObj = graph(adjMat);
L = distances(graphObj); % Calculate the shortest path length
%Average Path Length of the graph
APL = sum(sum(L))/(size(L,1)*size(L,1)-1);
end

在这两个脚本的基础上，我们可以对不同概率$p$下的网络结构进行属性统计，并绘制成图2。

%% 本脚本将被用于绘制小世界网络图的两个结构属性随边重连概率p变化的曲线图
% 尝试对文献“collective dynamics of small world networks”的图二复现
% 依赖的外部脚本：
% 1. WattsStrogatz，用于生成小世界网络
% 2. avgClusteringCoefficient，用于统计网络的群聚系数
% 3. avgPathLength，用于统计网络的平均路径长度
% 特别地，在可视化工作中，会用到第三方颜色库cbrewer
%% 代码实现
clear;
% 这个刻度需要好好调整
p = [0.00014,0.0003,0.00045,0.001,0.0019,0.004,0.008,0.018,0.021,0.071,0.13,0.23,0.5,1]; % 接近论文的概率组合
caseNum = length(p); % 仿真参数的组合个数
N = 1000;
k=5;
%% 计算regular ring的两个结构参数
regularRing = WattsStrogatz(N,k,0); % 生成得到regular ring
regularAM = adjacency(regularRing); % 获得邻接矩阵
C0 = avgClusteringCoefficient(regularAM);
APL0 = avgPathLength(regularAM);
%% 开始统计不同概率下得到的无序网络同规则网络的两个结构属性的比值
ratioC = zeros(caseNum,1);
ratioAPL = zeros(caseNum,1);
C = ratioC;
APL = ratioAPL;
for i=1:caseNum
    graphI = WattsStrogatz(N,k,p(i));
    graphAM = adjacency(graphI); % 邻接矩阵
    C(i,1) = avgClusteringCoefficient(graphAM);
    APL(i,1) = avgPathLength(graphAM);
    ratioC(i,1) = avgClusteringCoefficient(graphAM)/C0;
    ratioAPL(i,1) = avgPathLength(graphAM)/APL0;
end % 更新演化结束
%% 绘制半对数图
figure
semilogx(p,ratioC,'o');
hold on;
semilogx(p,ratioAPL,'*');
%% 绘制双刻度y的双曲线图
% 参考来源：阿昆的科研笔记
%% 图片尺寸设置（单位：厘米）
figureUnits = 'centimeters';
figureWidth = 20;
figureHeight = 8;
%% 定义线型、标记符号、线宽和颜色
colorTypes = cbrewer('div','Spectral',8,'linear');
% 定义因变量Y1线型、符号、线宽与颜色
Y1_LS = '-';
Y1_MK = 'none';
Y1_LW = 1.5;
Y1_C = colorTypes(2,:);
% 定义因变量Y2线型、符号、线宽与颜色
Y2_LS = 'o';
Y2_MK = 'none';
Y2_LW = 1.5;
Y2_C = colorTypes(7,:);
%% 窗口设置
figure
%% 激活左轴
yyaxis left
semilogx(p,C,'o','Color', Y1_C,'LineWidth',1.5);
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight],'Color','w'); 
hold on;
hYLabel1 = ylabel('');
% 坐标区调整
set(gca, 'YColor', [.1,.1,.1],...         % 坐标轴颜色
        'YTick', 0:0.2:1,...        % 刻度位置、间隔
        'Ylim' , [0 1], ...          % 坐标轴范围
        'Yticklabel',{[0:0.2:1]},...
        'XMinorGrid','on',...
        'YMinorGrid','on',...
        'YMinorTick','on');    % Y坐标轴刻度标签
ylabel('Cluster Coefficient');
%% 激活右轴
yyaxis right
semilogx(p,APL,'s','Color', Y2_C,'LineWidth',1.5)
hold on;
% 坐标区调整
set(gca, 'YColor', [.1,.1,.1],...         % 坐标轴颜色
        'YTick', 0:10:55,...        % 刻度位置、间隔
        'Ylim' , [0 55], ...          % 坐标轴范围
        'Yticklabel',{0:10:55},...
        'XMinorGrid','on',...
        'YMinorGrid','on',...
        'YMinorTick','on',...
        'MinorGridAlpha',0.4);    % Y坐标轴刻度标签

legend('Cluster Coefficient', 'Average Path Length');
xlabel('$p$','Interpreter','latex');
ylabel('Average Path Length');
xticklabels({0.0001,0.001,0.01,0.1,1});
 %% 图片输出
%  fileout = 'test';
%  print(figureHandle,[fileout,'.png'],'-r600','-dpng');

拓展补充

在绘制群聚系数和平均路径长度随边重连概率$p$的过程中，我们将这两个结构属性通过双轴散点图进行了直观地展示。
在这一过程中，可以拓展了解一下

1. matlab着色方案cbrewer。
2. 双轴散点图/曲线图绘制方案
   

参考资料/文献

[1]: D. J. Watts, S. H. Strogatz. Collective dynamics of ‘small­world’ networks[J]. nature, 1998, 393(6684): 440­442
[2]: Build Watts-Strogatz Small World Graph Model