一、导言
生成树（spanning tree）：在图论中，无向图G=(V,E)的生成树（spanning tree)是具有G的全部顶点，但边数最少的联通子图。假设G中一共有n个顶点，一颗生成树满足下列条件：
（1）n个顶点；
（2）n-1条边；
（3）n个顶点联通；
（4）一个图的生成树可能有多个。
最小生成树（minimum spanning tree， MST）/最小生成森林：联通加权无向图中边缘权重加和最小的生成树。给定无向图G=(V,E)，（u,v）代表顶点u与顶点v的边，w(u,v)代表此边的权重，若存在生成树T使得：
最小，则T为G的最小生成树。对于非连通无向图来说，它的每一连通分量同样有最小生成树，它们的并被称为最小生成森林。最小生成树除了继承生成树的性质之外，还存在下面两个特点：
- 当图的每一条边的权值都相同时，该图的所有生成树都是最小生成树；
- 如果图的每一条边的权值都互不相同，那么最小生成树将只有一个。
最小生成树示例：
图1. 加权无向联通图及其最小生成树
最小生成树MST的实际应用：
- 构建成本最小的连接网络：电网，计算机网络、交通网络、供应链网络；
-
最小生成树聚类：考虑数据集D，计算两点之间的距离当作边缘权重（一般采用欧式距离）。最小生成树聚类思想为，首先通过Prim等算法对数据集生成最小生成树，然后根据给定的阈值对最小生成树的边权重进行扫描，当边缘权重大于设定的阈值时，将对应的边切断。对所有数据执行上述操作后，剩下的还相互连接的部分即为相同的类或簇。
图2. 基于MST的聚类示例
- 反洗钱：核心交易网络发现+核心交易路线分
二、Prim算法介绍
Prim's Algorithm也翻译做普里姆算法，是1930年捷克数学家算法沃伊捷赫·亚尔尼克发现；并在1957年由美国计算机科学家罗伯特·普里姆独立发现；1959年，艾兹格·迪科斯彻再次发现了该算法。因此，在某些场合，普里姆算法又被称为DJP算法、亚尔尼克算法或普里姆－亚尔尼克算法。
Prim算法的思想：从任意一个顶点开始，每次选择与当前顶点最近的一个顶点，并将两点之间的边加入到树中。算法描述如下：
- 输入：加权无向图G， 顶点集合为V，边集合为E；
- 初始化:
，s为集合V中选择的任意起始点，
；
- 重复下列操作，直到
:
3.1 在集合E中选取权值最小的边（u, v），其中u为集合中的元素，而v则是V中没有加入
3.2 将v加入集合中；
- 输出：使用集合
下面对算法的图例描述：
输入：
原始加权无向连通图
说明： 此为原始加权无向连通图，每条边的数字代表其权重；
第1步：顶点A、B、E和F通过单条边与D相连。A是距离D最近的顶点，因此将A及对应边AD以高亮表示
第2步：下一个顶点为距离D或A最近的顶点。B距D为9，距A为7，E为15，F为6。因此，F距D或A最近，因此将顶点F与相应边DF以高亮表示。
第3步：算法继续重复上面的步骤。距离A为7的顶点B被高亮表示。
第4步：在当前情况下，可以在C、E与G间进行选择。C距B为8，E距B为7，G距F为11。E最近，因此将顶点E与相应边BE高亮表示。
第5步：这里，可供选择的顶点只有C和G。C距E为5，G距E为9，故选取C，并与边EC一同高亮表示。
第6步：顶点G是唯一剩下的顶点，它距F为11，距E为9，E最近，故高亮表示G及相应边EG。
输出：所有顶点均已被选取，图中绿色部分即为连通图的最小生成树。在此例中，最小生成树的权值之和为39。
三、Kruskal算法介绍
Kruskal是另一个计算最小生成树的算法，由Joseph Kruskal在1956年发表，属于贪婪算法。
算法原理：将每个顶点放入其自身的数据集中，然后按照权重的升序来选择边。当选择每条边时，判断定义边的顶点是否在不同的数据集中。如果是，将此边插入最小生成树的集合中，将集合中包含每个顶点的联合体取出。算法描述如下：
- 新建图G，G中拥有原图中相同的节点，但不包含边；
- 将原图中所有的边按权值从小到大排讯；
- 从权值最小的边开始，如果这条边连接的两个节点在G中不在同一个分量，则添加这条边到图G中；
重复3，直到图G中所有的节点在同一个联通分量中
Kruskal算法的动态图示
Kruskal算法的演示如下：
假设有一张图G，包含若干点和边。
第1步：将所有的边的长度排序，用排序的结果作为选择边的依据。资源排序，对局部最优的资源进行选择，排序完成后，率先选择了边AD。
第2步：在剩下的边中寻找，找到了CE。这里边的权重也是5
第3步：依次类推找到了DF，AB，BE。
第4步：下面继续选择， BC或者EF尽管现在长度为8的边是最小的未选择的边。但是现在他们已经连通了（对于BC可以通过CE,EB来连接，类似的EF可以通过EB,BA,AD,DF来接连）。所以不需要选择他们。类似的BD也已经连通了（这里上图的连通线用红色表示了）。最后就剩下EG和FG了。当然我们选择了EG。（等判断是否联通）
四、Prim和Kruskal算法的R及python实现
1. R 实现 Prim机Kruskal算法：igraph包mst函数, RBGL包的mstree.prim/mstree.kruskal
RBGL代码示例：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Rgraphviz")
BiocManager::install("RGBL")
require(RGBL)
require(XML)
require(Rgraphviz)
require(ape) # DNA sequences
# 1. simple DNA example
cat(">No305",
"NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACT",
">No304",
"ATTCGAAAAACACACCCACTACTAAAAATTATCAACCACT",
">No306",
"ATTCGAAAAACACACCCACTACTAAAAATTATCAATCACT",
">No307",
"ATTCGAATAACACAGCCACTTCTAAAAATTATCAATCACT",
file = "exdna.fas", sep = "\n")
data <- read.dna("exdna.fas", format = "fasta")
#generate a distance matrix
dist <- dist.dna(data,model="raw", as.matrix=TRUE)
#creates an undirected graph
dist.g<-as(dist,Class="graphNEL")
#generates the minimum spanning tree using the prim algorithm
ms<-mstree.prim(dist.g)
fromto<-cbind(ms$edgeList[1,],ms$edgeList[2,],ms$weight[1,])
adjMST<-ftM2adjM(as.matrix(fromto[,1:2]),W=fromto[,3],edgemode="undirected")
am.graph<-new("graphAM", adjMat=adjMST )
plot(am.graph, attrs = list(node = list(fillcolor = "lightblue"),edge = list(arrowsize=0.5)),"neato")
#2. more complicated DNA example
a <- "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/"
b <- "SRR062641.filt.fastq.gz"
URL <- paste0(a, b)
download.file(URL, b)
X <- read.fastq(b)
names(X)<-paste0("DNA_1",1:length(X))
dist <- dist.dna(X[1:30],model="raw", as.matrix=TRUE)
#creates an undirected graph
dist.g<-as(dist,Class="graphNEL")
#generates the minimum spanning tree using kruskal algorithm
ms<-mstree.kruskal(dist.g)
fromto<-cbind(ms$edgeList[1,],ms$edgeList[2,],ms$weight[1,])
adjMST<-ftM2adjM(as.matrix(fromto[,1:2]),W=fromto[,3],edgemode="undirected")
am.graph<-new("graphAM", adjMat=adjMST )
plot(am.graph, attrs = list(node = list(fillcolor = "lightblue"),edge = list(arrowsize=0.5)),"neato")
结果示例：
> ms
$edgeList
[,1] [,2] [,3] [,4]
from "No305" "No306" "No305" "No306"
to "No305" "No304" "No306" "No307"
$weights
[,1] [,2] [,3] [,4]
weight 0 0.02631579 0.02631579 0.07894737
$nodes
[1] "No305" "No304" "No306" "No307"
结果1.png
> ms
$edgeList
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
from "DNA_114" "DNA_118" "DNA_119" "DNA_18" "DNA_11" "DNA_15" "DNA_110"
to "DNA_128" "DNA_124" "DNA_118" "DNA_14" "DNA_129" "DNA_117" "DNA_127"
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
from "DNA_120" "DNA_121" "DNA_122" "DNA_114" "DNA_119" "DNA_19" "DNA_120"
to "DNA_122" "DNA_14" "DNA_124" "DNA_11" "DNA_121" "DNA_126" "DNA_128"
[,15] [,16] [,17] [,18] [,19] [,20] [,21]
from "DNA_17" "DNA_111" "DNA_116" "DNA_12" "DNA_15" "DNA_110" "DNA_110"
to "DNA_129" "DNA_129" "DNA_118" "DNA_17" "DNA_115" "DNA_11" "DNA_16"
[,22] [,23] [,24] [,25] [,26] [,27] [,28]
from "DNA_110" "DNA_111" "DNA_112" "DNA_125" "DNA_125" "DNA_121" "DNA_113"
to "DNA_126" "DNA_112" "DNA_15" "DNA_130" "DNA_127" "DNA_123" "DNA_19"
[,29]
from "DNA_11"
to "DNA_13"
$weights
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
weight 0.5833333 0.59375 0.59375 0.625 0.6354167 0.6458333 0.6458333 0.6458333
[,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
weight 0.6458333 0.6458333 0.6458333 0.65625 0.65625 0.65625 0.65625 0.6666667
[,17] [,18] [,19] [,20] [,21] [,22] [,23]
weight 0.6666667 0.6770833 0.6770833 0.6770833 0.6770833 0.6770833 0.6770833
[,24] [,25] [,26] [,27] [,28] [,29]
weight 0.6770833 0.6770833 0.6770833 0.6875 0.6875 0.6979167
$nodes
[1] "DNA_11" "DNA_12" "DNA_13" "DNA_14" "DNA_15" "DNA_16" "DNA_17"
[8] "DNA_18" "DNA_19" "DNA_110" "DNA_111" "DNA_112" "DNA_113" "DNA_114"
[15] "DNA_115" "DNA_116" "DNA_117" "DNA_118" "DNA_119" "DNA_120" "DNA_121"
[22] "DNA_122" "DNA_123" "DNA_124" "DNA_125" "DNA_126" "DNA_127" "DNA_128"
[29] "DNA_129" "DNA_130"
结果2.png
2. Python实现Prim和Kruskal算法
Prim算法示例：
import random
import time
def random_matrix_genetor(vex_num=10):
'''
随机图顶点矩阵生成器
输入：顶点个数，即矩阵维数
'''
data_matrix=[]
for i in range(vex_num):
one_list=[]
for j in range(vex_num):
one_list.append(random.randint(1, 100))
data_matrix.append(one_list)
return data_matrix
def prim(data_matrix):
'''
prim 算法
'''
vex_num=len(data_matrix)
prims=[]
weights=[]
flag_list=[False]*vex_num
node=0
for i in range(vex_num):
prims.append(0)
weights.append(0)
flag_list[node]=True
for i in range(vex_num):
weights[i]=data_matrix[node][i]
for i in range(vex_num-1):
min_value=float('inf')
for j in range(vex_num):
if weights[j]!=float('inf') and weights[j]<min_value and not flag_list[j]:
min_value=weights[j]
node=j
if node==0:
return
flag_list[node]=True
for m in range(vex_num):
if weights[m]>data_matrix[node][m] and not flag_list[m]:
weights[m]=data_matrix[node][m]
prims[m]=node
return weights, prims
def main_test_func(vex_num=10):
'''
主测试函数
'''
start_time=time.time()
data_matrix=random_matrix_genetor(vex_num)
weights, prims=prim(data_matrix)
print(weights)
print(prims)
end_time=time.time()
return end_time-start_time
if __name__=='__main__':
data_matrix=random_matrix_genetor(10)
weights, prims=prim(data_matrix)
print(weights)
print(prims)
time_list=[]
print('----------------------------10顶点测试-------------------------------------')
time10=main_test_func(10)
time_list.append(time10)
print('----------------------------50顶点测试-------------------------------------')
time50=main_test_func(50)
time_list.append(time50)
print('----------------------------100顶点测试-------------------------------------')
time100=main_test_func(100)
time_list.append(time100)
print('---------------------------------时间消耗对比--------------------------------')
for one_time in time_list:
print(one_time)
输出：
[26, 11, 12, 36, 13, 1, 4, 4, 6, 4]
[0, 7, 3, 0, 5, 7, 5, 2, 1, 7]
----------------------------10顶点测试-------------------------------------
[71, 4, 39, 8, 20, 17, 4, 13, 10, 1]
[0, 3, 5, 4, 0, 6, 4, 4, 0, 8]
----------------------------50顶点测试-------------------------------------
[60, 4, 4, 3, 4, 3, 9, 7, 1, 1, 3, 4, 2, 6, 5, 4, 5, 1, 2, 4, 3, 3, 2, 5, 4, 4, 4, 1, 4, 3, 5, 6, 2, 3, 5, 5, 3, 4, 2, 3, 4, 6, 3, 2, 2, 3, 4, 3, 3, 3]
[0, 27, 21, 37, 29, 1, 17, 16, 27, 27, 9, 21, 37, 23, 13, 37, 49, 34, 33, 26, 2, 42, 27, 40, 47, 16, 33, 15, 35, 26, 11, 40, 24, 24, 12, 13, 49, 11, 29, 19, 48, 7, 4, 7, 49, 5, 15, 0, 10, 45]
----------------------------100顶点测试-------------------------------------
[68, 1, 2, 2, 1, 3, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1, 3, 3, 1, 2, 3, 1, 2, 1, 1, 2, 1, 2, 3, 1, 1, 1, 2, 1, 2, 1, 2, 1, 1, 1, 3, 2, 2, 2, 3, 1, 1, 1, 2, 1, 3, 1, 2, 2, 1, 2, 2, 4, 2, 2, 1, 1, 1, 1, 3, 1, 2, 1, 3, 2, 2, 2, 1, 1, 1, 8, 5, 1, 2, 1, 2, 3, 3, 1, 2, 3, 6, 1, 1, 1, 1, 1, 2, 1, 1, 2, 2]
[0, 4, 7, 1, 42, 99, 97, 11, 68, 21, 28, 27, 50, 26, 49, 6, 7, 22, 44, 80, 95, 3, 34, 30, 41, 49, 35, 52, 31, 12, 0, 63, 22, 90, 59, 15, 86, 74, 36, 82, 27, 32, 91, 45, 2, 82, 91, 55, 66, 34, 59, 86, 63, 52, 82, 56, 0, 82, 98, 0, 7, 38, 59, 26, 39, 31, 37, 56, 65, 87, 57, 2, 62, 56, 94, 28, 46, 38, 37, 29, 68, 14, 38, 54, 54, 84, 14, 6, 42, 27, 57, 12, 51, 4, 29, 12, 95, 30, 47, 35]
---------------------------------时间消耗对比--------------------------------
0.0002651214599609375
0.004545927047729492
0.020392179489135742
Kruskal算法示例：
node = dict()
rank = dict()
def make_set(point):
node[point] = point
rank[point] = 0
def find(point):
if node[point] != point:
node[point] = find(node[point])
return node[point]
def merge(point1, point2):
root1 = find(point1)
root2 = find(point2)
if root1 != root2:
if rank[root1] > rank[root2]:
node[root2] = root1
else:
node[root1] = root2
if rank[root1] == rank[root2] : rank[root2] += 1
def Kruskal(graph):
for vertice in graph['vertices']:
make_set(vertice)
mst = set()
edges = list(graph['edges'])
edges.sort()
for edge in edges:
weight, v1, v2 = edge
if find(v1) != find(v2):
merge(v1 , v2)
mst.add(edge)
return mst
graph = {
'vertices': ['A', 'B', 'C', 'D','E'],
'edges': set([
(1, 'A', 'B'),
(5, 'A', 'C'),
(3, 'A', 'D'),
(4, 'B', 'C'),
(2, 'B', 'D'),
(1, 'C', 'D'),
(8, 'B', 'E'),
(3, 'A', 'E'),
])
}
print(Kruskal(graph))
输出：
{(1, 'C', 'D'), (1, 'A', 'B'), (3, 'A', 'E'), (2, 'B', 'D')}
参考资料
[1] Spanning tree, wiki: https://zh.wikipedia.org/wiki/%E6%9C%80%E5%B0%8F%E7%94%9F%E6%88%90%E6%A0%91;
[2]Minimum spanning tree, wiki: https://zh.wikipedia.org/wiki/%E7%94%9F%E6%88%90%E6%A0%91;
[3] 最小生成树聚类: http://dataminer.me/2017/10/20/%E6%9C%80%E5%B0%8F%E7%94%9F%E6%88%90%E6%A0%91%E8%81%9A%E7%B1%BB/;
[4] 最小生成树-Prim算法和Kruskal算法: https://www.cnblogs.com/biyeymyhjob/archive/2012/07/30/2615542.html;
[5] prim's algorithm, wiki: https://zh.wikipedia.org/wiki/%E6%99%AE%E6%9E%97%E5%A7%86%E7%AE%97%E6%B3%95；
[6]图的基本算法（最小生成树）: https://www.jianshu.com/p/efcd21494dff;
[7]算法导论--最小生成树（Kruskal和Prim算法）:https://blog.csdn.net/luoshixian099/article/details/51908175;
[8]Kruskal算法，wiki：https://zh.wikipedia.org/wiki/%E5%85%8B%E9%B2%81%E6%96%AF%E5%85%8B%E5%B0%94%E6%BC%94%E7%AE%97%E6%B3%95;
[9] RGBL installation: http://www.bioconductor.org/packages/release/bioc/html/RBGL.html；
[10]mstree.kruskal: https://www.rdocumentation.org/packages/RBGL/versions/1.48.1/topics/mstree.kruskal;
[11]Draw Minimum Spanning Tree (Mst) With Pie Charts In R: https://www.biostars.org/p/18876/;
[12]graphviz installation:https://www.bioconductor.org/packages/release/bioc/html/Rgraphviz.html;
[13]ape, read.dna 函数：https://www.rdocumentation.org/packages/ape/versions/5.3/topics/read.dna;
[14]Minimum spanning tree in igraph: https://igraph.org/r/doc/mst.html;
[15]python实现Prim算法求解加权连通图的最小生成树问题:https://blog.csdn.net/Together_CZ/article/details/74783631;
[16] Github: https://github.com/qiwsir/algorithm/blob/master/kruskal_algorithm.md;
[17]知乎：https://zhuanlan.zhihu.com/p/61628249
网友评论