Matlab实现小世界网络生成及其分析
- 引言
- 社会网络的结构属性
- Watts--Strogatz小世界网络生成算法
-
- matlab代码实现
- 网络结构分析
-
- 网络结构属性分析及实现代码
- 拓展补充
- 参考资料/文献
引言 本文的主要工作包括:
- 对社会网络的结构属性做了一个介绍
- 对群聚系数、平均路径长度以及节点度分布计算方式做了一个介绍
- 提供了Watts--Strogatz“小世界”网络生成算法流程及其matlab代码实现
进一步,为了描述不同网络结构的图的表示的特点,学者们提出并研究了许多图的属性特征。常见的图结构属性包括:
- 节点度:指与该节点相连的边的数目。在无向图中,节点度即为与该节点相连接的所有其他节点的数目。而在有向图中,将所有指向某节点的边的数量记为该节点的入度,所有从该节点出发指向别的节点的边的数量记为该节点的出度。网络平均度反应了网络的疏密程度。
- 度分布:指随机选定一个节点,该节点的度数等于 k k k的概率 p ( k ) p(k) p(k),即网络中节点的度数等于 k k k的比例。通过度分布则可以刻画不同节点的重要性。
- 群聚系数:用于描述网络中与同一节点相连的节点间也互为相邻节点的程度。其用于刻画社会网络中一个人朋友们之间也互相是朋友的概率,反映了社会网络中的聚集性。
- 介数:为图中某节点承载整个图所有最短路径的数量,通常用来评价节点的重要程度,比如在连接不同社群之间的中介节点的介数相对于其他节点来说会非常大,也体现了其在社会网络信息传递中的重要程度。
- 网络直径:指网络中任意两个节点之间的最短路径的最大值,影响着复杂系统的收敛时间和整个网络的稳定性。
- 平均路径长度:指网络中任意两个节点之间的最短路径的平均值。用来度量整体网络节点之间通信的有效性,即一个个体的信息平均传播多少距离才能被另一个个体接收。
Watts--Strogatz小世界网络生成算法 前面,我们已经了解到不同的社会网络会具有不同网络结构属性。并且,我们可以根据
群聚系数、平均路径长度
的差异将网络结构划分为三类:- 规则网络:具有明显的局部群聚结构和较大的网络群聚系数,但网络平均路径长度会随网络节点数线性增长。
- 随机网络:网络平均路径长度随网络节点数对数增长,但不具有局部群聚结构,即网络群聚系数很小。
- “小世界”网络:具有同规则网络相当的群聚系数以及和随机网络一样较小的网络平均路径长度。
基于此,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));
endretG = graph(s,t);
% 返回图
retMat = adjacency(retG);
% 返回邻接矩阵
delete(hwait);
end
网络结构分析 在这一重连过程中,网络结构会随着边重连概率 p p p的值的不同而有所改变。
文章图片
相应地,网络结构属性也会有所改变。我们只考虑网络结构的
群聚系数
以及平均路径长度
的改变情况。并且,统计了不同的边重连概率 p p p下的结构属性的变化情况。文章图片
图2:群聚系数与平均路径长度随边重连概率变化的曲线图
从图2中,我们可以发现:
- 当边重连概率 p p p在区间[0,0.01]增长时,网络结构的群聚系数不会发生太大的变化,而网络的
平均路径长度
却会快速下降。这就是“捷径边”的影响。影响了网络的全局结构属性 - 而当边重连概率进一步增加时,原来网络结构中的边随机性增强。最终在 p = 1 p=1 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 p 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 p p的过程中,我们将这两个结构属性通过双轴散点图
进行了直观地展示。在这一过程中,可以拓展了解一下
- matlab着色方案
cbrewer
。 - 双轴散点图/曲线图绘制方案
文章图片
[2]: Build Watts-Strogatz Small World Graph Model
推荐阅读
- 图像算法|非局部均值滤波算法(NL-means)
- dfs|dfs bfs连通区域算法 matlab,【算法】图论(一) —— 基本图算法(BFS/DFS/强连通分量)...
- 算法优化|回归预测 | MATLAB实现PSO-RBF多输入单输出
- MATLAB学习|MTALAB学习笔记——二三维图像的基本画法
- MATLAB学习|MATLAB运用——计算三维物体的质心(水花号)
- #|基于改进二进制粒子群算法的配电网重构
- 编程|基于Matlab的遗传算法在图像分割问题中的应用
- #|神奇的量子世界——量子遗传算法(Python&Matlab实现)
- 机器学习|遗传算法与matlab实现 Genetic Algorithm