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

山顶夕景

已于 2024-05-23 22:21:38 修改

阅读量7.5k

点赞数 14

分类专栏：数学math 文章标签： python 线性代数

于 2021-09-11 16:02:44 首次发布

本文链接：https://blog.csdn.net/qq_35812205/article/details/120234612

版权

数学math 专栏收录该内容

36 篇文章 44 订阅

订阅专栏

学习心得

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

文章目录

一、背景介绍

1.1 图论基础

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

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

一个无向无权图的例子：

其邻接矩阵为：
$\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)$

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

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

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

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

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

1.2 拉普拉斯矩阵的变体

节点数为 $n$ 的简单图 $G$ ， $A$ 是 $G$ 的邻接矩阵，则如上面介绍的那样， $G$ 的拉普拉斯矩阵即 $L = D - A$ 。
（1）对称归一化后的拉普拉斯矩阵： $L^{s y s}=D^{-1 / 2} L D^{-1 / 2}=I-D^{-1 / 2} A D^{-1 / 2}$
（2）随机游走归一化的拉普拉斯矩阵： $\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}$

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

拉普拉斯矩阵是半正定对称矩阵
对称矩阵有n个线性无关的特征向量，n是Graph中节点的个数 $\Rightarrow$ 拉普拉斯矩阵可以特征分解
半正定矩阵的特征值非负
对称矩阵的特征向量构成的矩阵为正交阵 $\Rightarrow U^{T} U=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-16j
 1.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变换

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

$f(x)=\int_{\mathbb{R}} \mathcal{F}(v) e^{2 \pi i x \cdot v} d v$

一阶导数定义
$f^{\prime}(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h}$
拉普拉斯相当于二阶导数
$\Delta f(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-2 f(x)+f(x-h)}{h^{2}}$
在graph上，定义一阶导数为
$f_{* g}^{\prime}(x)=f(x)-f(y)$
对应的拉普拉斯算子定义为
$\Delta_{*g} f^{\prime}(x)=\Sigma_{y \sim x} (f(x)-f(y))$
假设 $D$ 为 $N\times{N}$ 的度矩阵(degree matrix)
$j)=\left\{\begin{array}{ll}{d_{i}} & {\text { if } i=j} \\ {0} & {\text { otherwise }}\end{array}\right.$
$A$ 为 $N\times{N}$ 的邻接矩阵(adjacency matrix)
$j)=\left\{\begin{array}{ll}{1} & {\text { if } x_{i} \sim x_{j}} \\ {0} & {\text { otherwise }}\end{array}\right.$
那么图上的Laplacian算子可以写成
$L = D - A$
标准化后得到
$L=I_{N}-D^{-\frac{1}{2}} A D^{-\frac{1}{2}}$
定义Laplacian算子的目的是为了找到Fourier变换的基。

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

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

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

在这里插入图片描述

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