快捷导航

特征向量中心性计算中的迭代是怎么回事?

2021-9-8 11:53| 发布者: Fuller| 查看: 180| 评论: 0

摘要: 有同学在研读这篇论文《突发公共事件网络新闻文本集信息关联及可视化》,看到了特征向量中心性,网上搜到的大部分资料并没有仔细讲解计算公式是怎么来的,使用《Jupyter Notebook使用Python计算特征向量中心度(Eigen ...

有同学在研读这篇论文《突发公共事件网络新闻文本集信息关联及可视化》,看到了特征向量中心性,网上搜到的大部分资料并没有仔细讲解计算公式是怎么来的,使用《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

计算特征向量.zip


鲜花

握手

雷人

路过

鸡蛋

最新评论

GMT+8, 2021-9-22 16:27