有同学在研读这篇论文《突发公共事件网络新闻文本集信息关联及可视化》，看到了特征向量中心性，网上搜到的大部分资料并没有仔细讲解计算公式是怎么来的，使用《Jupyter Notebook使用Python计算特征向量中心度(Eigenvector Centrality)》介绍的networkx函数时，会问：为什么还要迭代？那么我们就从特征向量的计算开始探索，复习一下特征向量可以怎样计算。 使用Jupyter Notebook的原因是可以方便地记录下来整个探索过程，如果学习和整理知识的效率远比编程的效率重要，或者知识和数据探索的效率远比编程的效率重要，那么用Jupyter Notebook编程，同时记录探索和学习过程，是很适合的。 在这个Notebook中，我们用sympy程序包计算特征向量和矩阵运算，因为sympy是符号运算，跟numpy相比，通俗地理解就是sympy可以做精确的运算和推导，而numpy是数值运算，会有数制误差和截断误差累积下来。 1，求矩阵的特征向量 1.1，引入程序包 import sympy as sp sp.init_printing() 上面的init_printing可以把矩阵显示的更加好看。 参考Sympy官方手册，实验一下求特征向量，先从官方手册网页上拷贝一段代码。这是一个很“完美”的矩阵，实际应用场景几乎不可能遇到的，主要是为了演示原理。 1.2，用sympy函数计算特征向量 from sympy.matrices import Matrix M1 = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1]) M1
output1 = M1.eigenvects() output1
eigenvects()这个函数可以求出来特征向量和特征值，上面的结果中，第一个数是特征值，第二个数是乘数（multiplicity）或者重数，表示同一个特征值重复出现了几次。那么，还可以怎样计算特征向量呢？下面我们试一试 2，用networkx计算特征向量 2.1，构造图 使用M1矩阵构造一个图，有向边是用列向量描述的。 import networkx as nx import matplotlib.pyplot as plt import pylab import numpy as np vertices_s1=np.array([0, 0, 1, 1, 2, 2]) vertices_e1=np.array([1, 2, 0, 2, 0, 2]) value1=np.array([1, 1, 1, 1, 1, 1]) G1=nx.DiGraph() for i in range(np.size(vertices_s1)): G1.add_weighted_edges_from([(vertices_s1[i],vertices_e1[i],value1[i])]) 2.2，画图观察 pos=nx.shell_layout(G1) nx.draw(G1,pos,with_labels=True, node_color='white', edge_color='red', node_size=800, alpha=1 ) pylab.title('BetweennessCentralityTest',fontsize=25) pylab.show()
2.3，计算向量中心度 eigenvector1 = nx.eigenvector_centrality(G1, max_iter=5000) #特征向量中心度中心度 print("输出特征向量中心度的计算值：") for item in eigenvector1: print(item,"\t",eigenvector1[item]) 输出特征向量中心度的计算值： 0 0.5345226275102508 1 0.2672618166544383 2 0.8017834383660635 3，做一个矩阵乘向量的实验 随意构造一个初始向量吧，就是(1, 1, 1)。为了输入方便，在正文里这样写，约定为一个列向量。 X1 = Matrix([1, 1, 1]) X1
M1*X1
for i in range(10): X1 = M1 * X1 X1
n1 = X1.norm() X1 = X1 / n1 X1
normalize一个向量应该有专门的函数吧？但是我没有找到，所以就手工写代码用上面两步实现的。 实验一下再乘1000次以后会怎样 for i in range(1000): X1 = M1 * X1 X1 = X1 / X1.norm() X1
是不是很神奇，这个值再也不变了。 下面我们再回过头来看看前面用sympy函数求出来的特征向量 ev1 = output1[2] ev1
我们只看其中的特征向量 ev1[2][0]
执行normalization ev1_norm = ev1[2][0] / ev1[2][0].norm() ev1_norm
迭代乘出来的和用求特征向量函数算出来的竟然是一样的。 变成小数看看，方便后面的对比观察。发现跟networkx计算的结果也是一样。 print("x1=", float(ev1_norm[0]), "; x2=", float(ev1_norm[1]), "; x3=", float(ev1_norm[2])) 输出结果是： x1= 0.5345224838248488 ; x2= 0.2672612419124244 ; x3= 0.8017837257372732 4，换一个初始向量再做一次迭代乘实验 X2 = Matrix([1, 1, 0]) X2
for i in range(1000): X2 = M1 * X2 X2 = X2 / X2.norm() X2
上面这个运算好花时间，搞得电脑风扇嗡嗡响了好几分钟，很奇怪，这个为什么这么慢？可能是符号运算的原因，符号运算计算量太大了，等使用numpy做个比较就知道数值计算有多快了。 我都想现在中断这个运算，但是不知道怎么中断，因为我想改变一下运算顺序。百度上搜一下。有人说ctrl+C，试了不行。其实本页面最顶上有个方块图案的按钮，那就是中断。 改成下面的程序，把normalize函数放在循环外边，最后再执行，这样能省很多。 for i in range(1000): X2 = M1 * X2 X2 = X2 / X2.norm() X2
如果你真的执行了上面的代码，会看到一个好长的输出结果。符号运算太吊了，竟然算出一个那么大的分数来，怪不得这么慢，试试得到每个元素的值是什么 print("x1=", float(X2[0]), "; x2=", float(X2[1]), "; x3=", float(X2[2])) 输出结果是： x1= 0.5345224838248488 ; x2= 0.2672612419124244 ; x3= 0.8017837257372732 费了老大劲，还是第三个特征向量，找个什么样的初始向量能迭代出第一个特征向量呢？看来应该第一个元素是负的 X3 = Matrix([-3, 5, 1]) X3
for i in range(100): X3 = M1 * X3 X3 = X3 / X3.norm() X3
这次迭代次数少一些，别让风扇那么累，看起来还是差不多 print("x1=", float(X3[0]), "; x2=", float(X3[1]), "; x3=", float(X3[2])) 输出结果是： x1= 0.5345224838248488 ; x2= 0.2672612419124244 ; x3= 0.8017837257372732 第三个元素置0吧 X4 = Matrix([-3, 5, 0]) X4
for i in range(100): X4 = M1 * X4 X4 = X4 / X4.norm() X4
print("x1=", float(X4[0]), "; x2=", float(X4[1]), "; x3=", float(X4[2])) 输出结果是： x1= 0.5345224838248488 ; x2= 0.2672612419124244 ; x3= 0.8017837257372732 想搞出一个不同的特征向量来，搞不出来啊，干脆就用第一个特征向量来迭代，看看还出不来吗？ X5 = Matrix([-1, 1, 0]) for i in range(100): X5 = M1 * X5 X5 = X5 / X5.norm() X5
这下搞出一个不同的特征向量了，试试稍微改变一下 X6 = Matrix([-1, 2, 0]) for i in range(100): X6 = M1 * X6 X6 = X6 / X6.norm() X6
print("x1=", float(X6[0]), "; x2=", float(X6[1]), "; x3=", float(X6[2])) 输出结果是： x1= 0.5345224838248488 ; x2= 0.2672612419124244 ; x3= 0.8017837257372732 看来第三个特征向量太强了，就像一个大个的星球，吸引力太大了，经过迭代，大部分都让它吸引去了。要记住，从理论上说：由于特征值有大于1的，有小于1的，那么经过迭代产生的向量的方向都会逼近特征值大于1的特征向量上，除非初始向量就是跟小于1的特征值的特征向量同方向。 5，总结 David C. Lay, Steven R. Lay, Judi J. McDonald写的Linear Algebra and Its Applications这本书，第五章讲Eigenvalue和Eigenvector，在introduction部分讲了一个花斑猫头鹰生态的例子，这是一个动力系统（Dynamic System），随着时间的演进，本质上就是一次次迭代，最后应该稳定到一个状态。特征向量就是代表最后稳定到的状态。这一章甚至整本书都是特别好读的，适合我们这种搞工程的，现在有必要复习一下了： 1. 第一节Introductory Example 已经说过了，读起来挺有意思 2. 第5.6节 Discrete Dynamic System 解释了为什么会这样，核心是一个公式：
如果特征值有大于1的，还有小于1的，大于1的k次方就更大了，小于1的k次方就接近0了，就类似这个鞍形图的走势，最后都趋近于特征值最大的而且大于1的特征值对应的特征向量方向。当然，还有全部小于1的，还有全部大于1的，那就不是鞍形图了，要么都一直往外发散，要么都收缩到0 3. 第5.8节 Iterative Estimates for Eigenvalues 解释了怎样用迭代的方法求特征值和特征向量，其实只能求一个绝对占优的特征值和特征向量，就像这个例子，特征值是2，而且很明显比其他的大，那么试了那么多初始向量，都可以快速趋近特征向量。另外，这一章也讲了，如果没有一个绝对占优的怎么办。当然，为了计算出来特征值，一开始肯定不知道有没有占有的，这些复杂情况该怎么办。 那么，回到本notebook的题目：特征向量中心性计算中的迭代是怎么回事？是不是已经清楚了。回到《Jupyter Notebook使用Python计算特征向量中心度(Eigenvector Centrality)》这篇去看，第二个例子依然很神奇，那个矩阵的特征值是0，重数是8，显然很快就会迭代到0，为什么networkx还能算出来一个说的过去的结果？根据networkx的手册,eigenvector_centrality()函数用了Perron–Frobenius theorem，但是具体怎么实现的这个函数，似乎还要研究研究。 6，下载本notebook |
