Python求解拉普拉斯矩阵及其特征值

学习心得

（1）laplacian matrix就是无向图中定义 L=D−AL=D-AL=D−A，其中D为邻接矩阵，A为度矩阵（是一个对角矩阵）。本文用的python计算拉普拉斯矩阵及其特征值、特征向量。
（2）numpy中文文档：https://www.numpy.org.cn/user/，可以配spyder看对应变量的矩阵元素。
（3）求图拉普拉斯矩阵的特征值：时间复杂度是O(n3)O(n^3)O(n3)（n为节点数），所以当图很大时用该方法效率很低。

文章目录

学习心得
一、背景介绍
- 1.1 图论基础
- 1.2 拉普拉斯矩阵的变体
- 1.3 拉普拉斯矩阵的优良性质：
- 1.4 GCN为啥要用拉普拉斯矩阵
二、Python代码实现
三、关于图Fourier变换
Reference

一、背景介绍

1.1 图论基础

定义一（图的邻接矩阵）：

给定一个图G={V,E}\mathcal{G}=\{\mathcal{V}, \mathcal{E}\}G={V,E}，其对应的邻接矩阵被记为A∈{0,1}N×N\mathbf{A} \in\{0,1\}^{N \times N}A∈{0,1}N×N。Ai,j=1\mathbf{A}_{i, j}=1Ai,j=1表示存在从结点viv_ivi到vjv_jvj的边，反之表示不存在从结点viv_ivi到vjv_jvj的边。
在无向图中，从结点viv_ivi到vjv_jvj的边存在，意味着从结点vjv_jvj到viv_ivi的边也存在。因而无向图的邻接矩阵是对称的。
在无权图中，各条边的权重被认为是等价的，即认为各条边的权重为111。
对于有权图，其对应的邻接矩阵通常被记为W∈{0,1}N×N\mathbf{W} \in\{0,1\}^{N \times N}W∈{0,1}N×N，其中Wi,j=wij\mathbf{W}_{i, j}=w_{ij}Wi,j=wij表示从结点viv_ivi到vjv_jvj的边的权重。若边不存在时，边的权重为000。

一个无向无权图的例子：

其邻接矩阵为：
A=(0101110100010011000110110)\mathbf{A}=\left(\begin{array}{lllll} 0 & 1 & 0 & 1 & 1 \\ 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 \end{array}\right) A=⎝⎜⎜⎜⎜⎛0101110100010011000110110⎠⎟⎟⎟⎟⎞

定义二（拉普拉斯矩阵，Laplacian Matrix）：

给定一个图G={V,E}\mathcal{G}=\{\mathcal{V}, \mathcal{E}\}G={V,E}，其邻接矩阵为AAA，其拉普拉斯矩阵定义为L=D−A\mathbf{L=D-A}L=D−A，其中度矩阵 D=diag(d(v1),⋯,d(vN))\mathbf{D=diag(d(v_1), \cdots, d(v_N))}D=diag(d(v1),⋯,d(vN))。

定义三（对称归一化的拉普拉斯矩阵，Symmetric normalized Laplacian）：

给定一个图G={V,E}\mathcal{G}=\{\mathcal{V}, \mathcal{E}\}G={V,E}，其邻接矩阵为AAA，其规范化的拉普拉斯矩阵定义为

L=D−12(D−A)D−12=I−D−12AD−12\mathbf{L=D^{-\frac{1}{2}}(D-A)D^{-\frac{1}{2}}=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}} L=D−21(D−A)D−21=I−D−21AD−21

1.2 拉普拉斯矩阵的变体

节点数为nnn的简单图GGG，AAA是GGG的邻接矩阵，则如上面介绍的那样，GGG的拉普拉斯矩阵即L=D−AL=D-AL=D−A。
（1）对称归一化后的拉普拉斯矩阵：Lsys=D−1/2LD−1/2=I−D−1/2AD−1/2L^{s y s}=D^{-1 / 2} L D^{-1 / 2}=I-D^{-1 / 2} A D^{-1 / 2} Lsys=D−1/2LD−1/2=I−D−1/2AD−1/2
（2）随机游走归一化的拉普拉斯矩阵：Lrw=D−1L=I−D−1ALi,jrw:={1i=j−1diag⁡(vi)if i≠jand viadjacent to vj0othewise \begin{aligned} &L^{r w}=D^{-1} L=I-D^{-1} A \\ &L_{i, j}^{r w}:= \begin{cases}1 & \mathrm{i}=\mathrm{j} \\ \frac{-1}{\sqrt{\operatorname{diag}\left(v_{i}\right)}} & \text { if } i \neq j \text { and } v_{i} \text { adjacent to } v_{j} \\ 0 & \text { othewise }\end{cases} \end{aligned} Lrw=D−1L=I−D−1ALi,jrw:=⎩⎪⎪⎨⎪⎪⎧1diag(vi)−10i=j if i=j and vi adjacent to vj othewise

1.3 拉普拉斯矩阵的优良性质：

拉普拉斯矩阵是半正定对称矩阵
对称矩阵有n个线性无关的特征向量，n是Graph中节点的个数 ⇒\Rightarrow⇒ 拉普拉斯矩阵可以特征分解
半正定矩阵的特征值非负
对称矩阵的特征向量构成的矩阵为正交阵 ⇒UTU=E\Rightarrow U^{T} U=E⇒UTU=E

1.4 GCN为啥要用拉普拉斯矩阵

拉普拉斯矩阵可以谱分解（特征分解）GCN是从谱域的角度提取拓扑图的空间特征的。
拉普拉斯矩阵只在中心元素和一阶相邻元素处有非零元素，其他位置皆为0.
传统傅里叶变换公式中的基函数是拉普拉斯算子，借助拉普拉斯矩阵，通过类比可以推导出Graph上的傅里叶变换公式。

二、Python代码实现

这是networkX库对稀疏矩阵A的处理方式，有少量改进（将所有内容保持稀疏）：

n,m = A.shape
diags = A.sum(axis=1)
D = sps.spdiags(diags.flatten(), [0], m, n, format='csr')
D - A

numpy和scipy两个库中模块中都提供了线性代数的库linalg，但scipy的linalg的更全面一些。

# laplacian矩阵
import numpy as np
def unnormalized_laplacian(adj_matrix):# 先求度矩阵R = np.sum(adj_matrix, axis=1)degreeMatrix = np.diag(R)return degreeMatrix - adj_matrix# 对称归一化的laplacian矩阵
def normalized_laplacian(adj_matrix):R = np.sum(adj_matrix, axis=1)R_sqrt = 1/np.sqrt(R)D_sqrt = np.diag(R_sqrt)I = np.eye(adj_matrix.shape[0])return I - np.matmul(np.matmul(D_sqrt, adj_matrix), D_sqrt)

matlab的代码实现为：

L = diag(sum(A,2)) - A   % or L=diag(sum(A))-A because A is symmetric

那么如何求矩阵的特征值呢——利用numpy的linalg.eig：

# !/usr/bin/python
## -*-coding:utf-8 -*-import numpy as np
A = np.array([[3,-1],[-1,3]])
print('打印A：\n{}'.format(A))
a, b = np.linalg.eig(A)
print('打印特征值a：\n{}'.format(a))
print('打印特征向量b：\n{}'.format(b))

得到特征值和特征向量：

打印A：
[[ 3 -1] [-1  3]]
打印特征值a：
[4. 2.]
打印特征向量b：
[[ 0.70710678  0.70710678] [-0.70710678  0.70710678]]

三阶矩阵试试，回归考研线代的题目：

import numpy as np
A = np.array([[2,-3,1],[1,-2,1],[1,-3,2]])
print('打印A：\n{}'.format(A))
a, b = np.linalg.eig(A)
print('打印特征值a：\n{}'.format(a))
print('打印特征向量b：\n{}'.format(b))

结果如下，看结果的特征向量矩阵，每一列代表一个一个特征向量，都是

打印A：
[[ 2 -3  1][ 1 -2  1][ 1 -3  2]]
打印特征值a：
[2.09037533e-15+0.00000000e+00j 1.00000000e+00+5.87474805e-16j1.00000000e+00-5.87474805e-16j]
打印特征向量b：
[[0.57735027+0.j         0.84946664+0.j         0.84946664-0.j        ][0.57735027+0.j         0.34188085-0.11423045j 0.34188085+0.11423045j][0.57735027+0.j         0.17617591-0.34269135j 0.17617591+0.34269135j]]

三、关于图Fourier变换

根据卷积原理，卷积公式可以写成
f∗g=F−1{F{f}⋅F{g}}f * g=\mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\} f∗g=F−1{F{f}⋅F{g}}
正、逆Fourier变换
F(v)=∫Rf(x)e−2πix⋅vdx\mathcal{F}(v)=\int_{\mathrm{R}} f(x) e^{-2 \pi i x \cdot v} d x F(v)=∫Rf(x)e−2πix⋅vdx

f(x)=∫RF(v)e2πix⋅vdvf(x)=\int_{\mathbb{R}} \mathcal{F}(v) e^{2 \pi i x \cdot v} d v f(x)=∫RF(v)e2πix⋅vdv

一阶导数定义
f′(x)=lim⁡h→0f(x+h)−f(x)hf^{\prime}(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h} f′(x)=h→0limhf(x+h)−f(x)
拉普拉斯相当于二阶导数
Δf(x)=lim⁡h→0f(x+h)−2f(x)+f(x−h)h2\Delta f(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-2 f(x)+f(x-h)}{h^{2}} Δf(x)=h→0limh2f(x+h)−2f(x)+f(x−h)
在graph上，定义一阶导数为
f∗g′(x)=f(x)−f(y)f_{* g}^{\prime}(x)=f(x)-f(y) f∗g′(x)=f(x)−f(y)
对应的拉普拉斯算子定义为
Δ∗gf′(x)=Σy∼x(f(x)−f(y))\Delta_{*g} f^{\prime}(x)=\Sigma_{y \sim x} (f(x)-f(y)) Δ∗gf′(x)=Σy∼x(f(x)−f(y))
假设DDD为N×NN\times{N}N×N的度矩阵(degree matrix)
D(i,j)={diif i=j0otherwise D(i, j)=\left\{\begin{array}{ll}{d_{i}} & {\text { if } i=j} \\ {0} & {\text { otherwise }}\end{array}\right. D(i,j)={di0 if i=j otherwise
AAA为N×NN\times{N}N×N的邻接矩阵(adjacency matrix)
A(i,j)={1if xi∼xj0otherwise A(i, j)=\left\{\begin{array}{ll}{1} & {\text { if } x_{i} \sim x_{j}} \\ {0} & {\text { otherwise }}\end{array}\right. A(i,j)={10 if xi∼xj otherwise
那么图上的Laplacian算子可以写成
L=D−AL=D-A L=D−A
标准化后得到
L=IN−D−12AD−12L=I_{N}-D^{-\frac{1}{2}} A D^{-\frac{1}{2}} L=IN−D−21AD−21
定义Laplacian算子的目的是为了找到Fourier变换的基。

传统Fourier变换的基就是Laplacian算子的一组特征向量
Δe2πix⋅v=λe2πix⋅v\Delta e^{2 \pi i x \cdot v}=\lambda e^{2 \pi i x \cdot v} Δe2πix⋅v=λe2πix⋅v
类似的，在graph上，有
Δf=(D−A)f=Lf\Delta{f}=(D-A)f=Lf Δf=(D−A)f=Lf
图拉普拉斯算子作用在由图节点信息构成的向量fff上得到的结果等于图拉普拉斯矩阵和向量fff的点积。

那么graph上的Fourier基就是LLL矩阵的nnn个特征向量U=[u1…un]U=\left[u_{1} \dots u_{n}\right]U=[u1…un]，LLL可以分解成L=UΛUTL=U \Lambda U^{T}L=UΛUT，其中，Λ\LambdaΛ是特征值组成的对角矩阵。

传统的Fourier变换与graph的Fourier变换区别

将f(i)f(i)f(i)看成是第iii个点上的signal，用向量x=(f(1)…f(n))∈Rnx=(f(1) \ldots f(n)) \in \mathbb{R}^{n}x=(f(1)…f(n))∈Rn来表示。矩阵形式的graph的Fourier变换为
GF{x}=UTx\mathcal{G F}\{x\}=U^{T} x GF{x}=UTx
类似的graph上的Fourier逆变换为
IGF{x}=Ux\mathcal{I} \mathcal{G} \mathcal{F}\{x\}=U x IGF{x}=Ux

Reference

（1）Python计算特征值与特征向量案例
（2）numpy的linalg官方文档：https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html
（3）用python求解特征向量和拉普拉斯矩阵
（4）图卷积网络(GCN)原理解析