谱图理论

写在前面

图神经网络（GNN）是近年来愈发火热，用以对图（Graph）数据进行特征提取以及各种下游任务。注意这里的图（Graph）应和图像（Image）区分，图是一种由点集与边集组成的数据结构，常记作\(\mathcal{G}=(\mathcal{V,E})\)。以下是我学习谱图理论时的一些记录。本文先介绍了线性代数的一些前置知识，包括内积与基、特征值与特征向量、二次型等内容。接着从基变换、频域信号的角度阐述了傅里叶变换的原理。然后通过分析拉普拉斯矩阵的特征值与特征向量，利用二次型的形式说明了其特征值及其对应的特征向量，和傅里叶变换中的频域信号之间的联系。最后使用图的数据结构表示图像，可视化其特征向量，直观感受图的低频和高频分量的差异。本人才疏学浅，错误难免，欢迎交流指正。

线性代数基础

内积与基

我们首先回顾一下线性代数的知识。本科第一次线性代数的时候，都是以繁杂的计算与证明为主，未曾有更直观、直觉的方式去理解。这里我想从线性变换的角度，讲述线性代数中与本文内容相关的知识。

我们定义实数向量 \(\vec x=[x_1, x_2,...,x_n]^T\) 与实数向量的 \(\vec x=[x_1, x_2,...,x_n]^T\) 内积为\(<\vec x,\vec y>=\sum x_iy_i\)。这是高中就学过的知识。但这里，我需要对向量的内积做一些扩充，即复数向量内积的定义。根据 hermitian 内积的定义^[1]，复平面内积定义为 \(<\vec x,\vec y>=\sum x_i\bar y_i\)。

如果我们将向量内积扩展到函数空间，我们可以将函数视为无限长的向量，如 \(f(x)=[f(x_1),f(x_2),...]\)，可以定义函数的内积^[2] \(<f,g>=\int_a^b f(x)g(x)dx\) 。若函数在复平面上，则可定义为 \(<f,g>=\int_a^b f(x)\overline{g(x)}dx\)（有时共轭会放在左边）。

然后我们来回顾一下向量空间中的基。由一组线性无关的向量可以张成一个向量空间，空间中的任意向量都可以使用这一组向量的线性组合表示。这组向量称作基向量。如果基向量构成的矩阵 \(A=[\vec\alpha_1,\vec\alpha_2,...,\vec\alpha_n]^T\) 满足 \(AA^T=E\)，则这组基向量称为正交基，若还满足 \(||\vec\alpha_i||=1,(i=1,2,...,n)\)，则称为标准正交基。如，二维向量 \([1,0]^T,[0,1]^T\) 构成一组标准正交基。在复平面中^[6]，转置被描述为 \((A^H)_{ij}=\overline{A_{ji}}\)，相比实数域，处了转置还要做一次共轭，称为共轭转置。那么正交基构成的矩阵 \(A\) 应满足 \(AA^H=E\) 。

假如有标准正交基\(\vec e_1, \vec e_2,...,\vec e_n\)，若该空间中向量 \[ \vec v=w_1\vec e_1+w_2\vec e_2+,...+w_n\vec e_n=[\vec e_1,\vec e_2,..., \vec e_n][w_1,w_2,...,w_m]^T \] 则 \[ [w_1,w_2,...,w_m]^T=[\vec e_1,\vec e_2,... ,\vec e_n]^{-1}\vec v=[\vec e_1,\vec e_2,... ,\vec e_n]^{T}\vec v \] 可以看到，若想获得一组某向量在一组标准正交基上的表示，只需要让该向量与这组正交基做内积即可。

函数也存在着基的概念，若函数 \(f(x)=a_ng(x)+b_nh(x)\)。那么 \(g(x)\) 和 \(h(x)\) 就是 \(f(x)\) 的基。若内积\(<g,h>=0\)，\(g\) 和 \(h\) 是一组正交基。

特征值与特征向量

简单地回顾一些特征值与特征向量的定义，设有矩阵\(A\)，若存在\(\lambda,\vec x\)，使得\(A\vec x=\lambda \vec x\)，则\(\lambda\)称为\(A\)的特征值，\(\vec x\)则称为特征向量，\(\lambda\)的值可以不止一个，同一个\(\lambda\)可以对应多个线性无关的\(\vec x\)。如果想了解特征值与特征向量的几何解释，可以参考^[5]，这不是本文的重点。

接下啦我想说明的是一种特殊的矩阵：实对称矩阵，即由实数组成的对称矩阵（若\(A=A^T\)则称\(A\)为对称矩阵。该矩阵有很多优良的性质。首先，\(N\) 阶是对称矩阵，具有 \(N\) 个特征值以及 \(N\) 个正交的特征向量。且实对称矩阵一定可以正交对角化（事实上，一个方阵能够正交对角化的充要条件就是该矩阵为对称矩阵），即 \[ Q^{-1}AQ=Q^TAQ=\Lambda \] 对角矩阵 \(\Lambda\) 对角线上为特征值 \(\lambda_1,\lambda_2,...\lambda_{N}\) ，分别对应 \(Q=[\vec q_1, \vec q_2,...,\vec{q}_{N}]\) 中的特征向量 \(\vec q_1,\vec q_2,...\vec{q}_{N}\) 。

由于其\(Q\)由一组线性无关的正交特征向量组成，他们可以构成一组正交基，特征向量组成的基也成为特征基。我们举个相似对角化的简单应用。假如要求实对称矩阵\(A\)的幂次预算，如\(A^{100}\)。直接计算是很麻烦的，若首先将其转化为对角矩阵 \(Q^{-1}AQ=\Lambda\)，那么 \(Q^{-1}AQ...Q^{-1}AQQ^{-1}AQ=Q^{-1}A^{100}Q=\Lambda^{100}\) 。则\(A^{100}=Q\Lambda^{100}Q^{-1}\) 。对角矩阵的100次幂是非常容易计算的，这就使得计算量大大减小。

二次型

二次型是线性代数中的一个重要概念，它涉及到二次多项式的表达和研究。在线性代数中，一个定义在 n 维向量空间上的二次型可以表示为一个二次多项式。一般地，一个n维向量 \(\mathbf{x} = [x_1, x_2, \ldots, x_n]^T\) 的二次型可表示为： \[ Q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x} \] 其中，\(A\) 是一个对称矩阵。对称矩阵的重要性在于它保证了二次型中交叉项的系数相等。在这个表达中，\(A\)是二次型的系数矩阵，\(\mathbf{x}\)是变量向量，\(\mathbf{x}^T\)表示向量的转置。

对于一个实二次型 \(Q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x}\)，其中 \(A\) 是对称矩阵，展开的形式如下： \[ Q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x} = \sum_{i=1}^{n} \sum_{j=1}^{n} a_{ij} x_i x_j \] 其中，\(a_{ij}\) 是矩阵 \(A\) 的元素。这个展开形式表示了二次型中所有可能的平方项和交叉项。

如果\(A\)是实对称矩阵，那么二次型称为实二次型。如果\(A\)是复对称矩阵，那么二次型称为复二次型。实二次型在应用中更为常见，因此下文主要讨论实二次型。

二次型 \(Q(\mathbf{x})\) 被称为正定的，如果对于所有非零的\(\mathbf{x}\)，都有\(Q(\mathbf{x}) > 0\)，如果在 0 处可以取等号，则称之为半正定的。负定的定义类似，只是符号相反。

实际应用中，二次型在优化问题、统计学、物理学等领域发挥着重要作用。例如，正定二次型常常出现在凸优化问题中，而负定二次型则在稳定性分析中有所应用。

注：该小结内容由 chatgpt 辅助编写。

谱分析

傅里叶变换

在学习傅里叶变换时，想必类似上面的图大家已经见过很多次了，大多数教材都会先从傅里叶级数入手，周期函数可以表示为许多正弦信号的叠加，有如下形式 \[ f(t)=\frac{a_0}{2}+\sum_{n=1}^{+\infty} [a_n\sin(\frac{n\pi}{l} t )+b_n\cos({\frac{n\pi}{l} t})] \] 其中，\(n\)是整数，\(l\)为半周期。\(1,\sin(n\omega t ),\cos({n\omega t})\) 可以看作是 \(f(t)\) 许许多多相互正交的基。这是因为^[3]：\(\int^{\pi}_{-\pi}\cos nxdx=0\),\(\int^{\pi}_{-\pi}\sin nxdx=0\), \(\int^{\pi}_{-\pi}\cos n_1x\sin n_2xdx=0\), \(\int^{\pi}_{-\pi}\cos n_1x\cos n_2xdx=0\),\(\int^{\pi}_{-\pi}\sin n_1x\sin n_2xdx=0\)。

因此我们可以将傅里叶变换看作是将时序信号使用不同频率的信号来表示，而这些信号之间是正交的如果使用 Euler 公式替换，可以转换为如下形式 \[ f(t)=\sum^{+\infty}_{n=-\infty}c_ne^{i\omega t} \] 其中，\(c_n=\frac{1}{l}\int_{-l}^{l}f(t)e^{-i\omega t}dt\)。\(l\)为半周期，\(\omega=\frac{n\pi}{l}\)，即频率（有些地方会将\(\omega\)视作角速度，指数项写作\(e^{2i\pi\omega t}\)）。这里 \(e^{i\omega t}\) 也是许多的标准正交基，这是因为 \(\int e^{i\omega_1 t}e^{-i\omega_2 t}=0\)，\(||e^{i\omega}||=1\)。那么当 \(l\rightarrow +\infty\) 时，\(c_n\) 就是傅里叶变换的形式了。即 \[ \mathcal{F}[f(t)]=\hat f(\omega)=\int_{-\infty}^{+\infty}f(t)e^{-i\omega t}dt \] 所以像函数 \(\hat f(\omega)\) 其实就是\(\omega\)对应的傅里叶基 \(e^{i\omega t}\) 的系数。上述式子也可以看作是 \(f(t)\) 与正交基 \(e^{i\omega t}\) 的内积，在内积与基中提到了，若某个向量希望使用某组正交基来表示，拿就让这个向量和这组标准正交基做内积（函数的内积公式见小节内积与基），就可以把基方向的分量提取出来了（仿佛一个筛子）。

拉普拉斯矩阵

我们假设存在图 \(\mathcal{G}=(\mathcal{V}, \mathcal{E})\)，结点 \(v\in \mathcal{V}\) 的特征为 \(x_{v}\in\mathcal{X}\)。为了简单起见，我们假设图为无向图，并且结点特征 仅有 1 个维度，因此所有的结点特征可以构成一个向量 \(x\)。

对于如上的图数据，它的拉普拉斯矩阵如下。其计算方法是\(D-A\)。\(D\)是每个结点的度组成的对角矩阵，\(A\)是图的邻接矩阵。 \[ L=\left(\begin{array}{rrrrrr} 2 & -1 & 0 & 0 & -1 & 0\\ -1 & 3 & -1 & 0 & -1 & 0\\ 0 & -1 & 2 & -1 & 0 & 0\\ 0 & 0 & -1 & 3 & -1 & -1\\ -1 & -1 & 0 & -1 & 3 & 0\\ 0 & 0 & 0 & -1 & 0 & 1\\ \end{array}\right) \] 由于 \(L\in \mathbb{R}^{n \times n}\) 是实对称、半正定矩阵，上一节中提到实对称矩阵的一些优良性质，存在标准正交基构成的矩阵 \(U\) 及其对应的特征值对角矩阵 \(\Lambda\) 使得拉普拉斯矩阵 \(L\) 满足 \[ U^TLU=\Lambda \] 由于拉普拉斯矩阵 \(L\) 为半正定矩阵（利用瑞丽熵即可证明），其对应的特征值按从小到大排序分别为 \(0=\lambda_{1},\lambda_{2},\dots,\lambda_{n}\) ，并构成对角矩阵 \(\Lambda\) ；其对应的特征向量构成矩阵 \(U=[ u_1, u_2,..., u_n]\) ，其中 \(\ u_i\) 为 \(n\) 个标准正交的特征向量。

我们暂且将 \(U\) 矩阵中的标准正交基看作是一组傅里叶基，那么我们希望得到各个基的分量，在内积与基以及傅里叶变换这两个小节都已提到相关的操作，类似的，就有以下等式： \[ \begin{gathered} x=U\hat{x}=\hat{x}_{1}u_{1}+\hat{x}_{2}u_{2}+\dots+\hat{x}_{n}u_{n} \\ \Rightarrow~ \hat{x}=U^Tx \end{gathered} \] \(\hat{x}\) 即为各个基的分量，或者用傅里叶变换的角度来说就是各个频域的分量。

问题来了，我们前面是暂且将它视为傅里叶基，那凭什么可以这么类比呢？我们来回忆一下傅里叶变换，它是将时序信号也是通过基变换，分解成了各个频率相互正交的频域信号的加权求和，而所谓的傅里叶基就是不同频率的频域信号。问题就在于，拉普拉斯矩阵 \(L\) 的特征向量，何以表示频域信号。

我们只需利用简单的变换，即可说明。 \[ u_{i}^TLu_{i}=u_{i}^T\lambda_{i} u_{i}=\lambda_{i} \] 其中，\(u_{i}^TLu_{i}\) 为二次型的形式。它可以表示什么呢？

根据小节二次型的公式，我们有以下公式推导： \[ \begin{align*} x^TLx&=\sum_{i\in \mathcal{V}}d_{i}x_{i}^2-\sum_{(i,j)\in \mathcal{E}}2x_{i}x_{j}\\ &=\sum_{(i,j)\in \mathcal{E}}(x_{i}^2-2x_{i}x_{j}+x_{j}^2)\\ &=\sum_{(i,j)\in \mathcal{E}}(x_{i}-x_{j})^2\\ \end{align*} \] 因此，二次型 \(x^TLx\) 可以衡量结点特征之间的差异。其结点之间差异越大，其震荡就越大，即越粗糙；否则，震荡越小，也就越平滑。

也就是说特征值 \(\lambda_{i}\) 可以体现特征向量的平滑程度。特征向量所对应的特征值越大，就说明震荡越明显，越不平滑，对应的是傅里叶变换中的高频信号；反之，则越平滑，对应傅里叶变换中的低频信号。

因此，我们将拉普拉斯矩阵的特征基视作一组傅里叶基是合理，我们利用拉普拉斯的特征分解，达到了和时序信号傅里叶变换类似的效果。我们将这种分解方法称为谱分解（Spectral Decomposition），相关的分析方式称为图的谱分析（Graph Spectral Analysis），相关理论称为谱图理论（Graph Spectral Theory）。

可视化

如何可视化呢图的频域信号呢？事实上，我们可以借助图像（Image）的结构来可视化。图像通常由一个矩阵表示，矩阵中的每一个值为像素点的值。根据通道数，图像又可以分为彩色图像和灰度图像。

事实上图像它也可表示称图的结构。如下图所示^[11]：

其中，每个像素与上下左右对角线相邻接。类似的，一个序列也可以转变为图的结构，同样也可以借助它来做可视化 。这里暂且利用图像的结构来可视化。

可视化的算法可描述如下：

• step1. 选择(10, 10)尺寸的图像结构，
• step2. 根据每个像素与上下左右对角线相邻接的规则构造邻接矩阵A, shape: (100, 100)
• step3. 根据邻接矩阵A构造拉普拉斯矩阵L, shape: (100, 100)
• step4. 对拉普拉斯矩阵进行对角化, 求得特征基矩阵U, shape(100, 100)，对角矩阵P, shape(100, 100)，注意上述2个矩阵均按特征值的大小升序排序。
• step5. 取特征值大小前10的特征基，每个特征基重构为图像结构Image, shape(10, 10)，进行可视化

经过以上的操作后，最后的结果如下图所示，特征值按从小到大排序。这就是每个频率所对应的分量。可以看到图像一开始是纯色的，后来开始出现了类似波形的渐变过程，随着频率的增加，其波形也越复杂。

以下为实现代码

import numpy as np
import math
import matplotlib.pyplot as plt

def getLaplacianOfImage(M, N):
    idx = np.array([i for i in range(M * N)]).reshape((M, N))
    tmp = np.ones((M+2, N+2), dtype=np.int32) * M * N
    tmp[1:M+1, 1:N+1] = idx

    directH = [1, -1, 0, 0, 1, -1, 1, -1]
    directV = [0, 0, 1, -1, 1, -1, -1, 1]
    A = np.zeros((M * N + 1, M * N + 1))
    for i in range(1, M+1):
        for j in range(1, N+1):
            for h, v in zip(directH, directV):
                A[tmp[i, j], tmp[i+h, j+v]] = 1
    A = A[0:M*N, 0:M*N]
    L = np.diag(np.sum(A, axis=1)) - A
    return L

M, N = 10, 10
L = getLaplacianOfImage(M, N)
eigenvalue, featurevector = np.linalg.eigh(L)
eigenbase = featurevector.T

min_val, max_val = np.min(eigenbase), np.max(eigenbase)
plt.figure(figsize=(15, 15))
for i in range(9):
    fv = eigenbase[i].reshape((M, N))
    plt.subplot(3, 3, i+1)
    plt.imshow(fv, vmin=min_val, vmax=max_val)
    plt.colorbar(fraction=0.045)

总结

写该博客的起因是学习GNN时，不理解为什么使用\(L\)矩阵，而不是用邻接矩阵或者其他实对称矩阵，以及在文献中^[13]看到 \(L\) 矩阵的谱分解其实就是傅里叶变换在图领域的应用。为了深入了解这二者的关系，捡起了以前学过但又未深入理解的知识。完成这篇博客，还是很有收获的。

参考资料

1. https://mathworld.wolfram.com/HermitianInnerProduct.html ↩︎
2. https://mathworld.wolfram.com/HilbertSpace.html ↩︎
3. 同济大学数学系．高等数学 [M]，第七版，高等教育出版社，2014-07-04 ↩︎
4. https://www.bilibili.com/video/BV1za411F76U?spm_id_from=333.337.search-card.all.click&vd_source=3eafcac5a31e0009a6433cea9bc7ab45 ↩︎
5. https://www.bilibili.com/video/BV1ys411472E?p=14&vd_source=3eafcac5a31e0009a6433cea9bc7ab45 ↩︎
6. https://math.mit.edu/~gs/linearalgebra/ ↩︎
7. https://en.wikipedia.org/wiki/Hilbert_space#Fourier_analysis ↩︎
8. https://en.wikipedia.org/wiki/Discrete_Fourier_transform ↩︎
9. https://en.wikipedia.org/wiki/DFT_matrix ↩︎
10. https://csustan.csustan.edu/~tom/Clustering/GraphLaplacian-tutorial.pdf ↩︎
11. https://distill.pub/2021/gnn-intro/ ↩︎
12. https://distill.pub/2021/understanding-gnns/#learning ↩︎
13. Li X, Zhu R, Cheng Y, et al. Finding Global Homophily in Graph Neural Networks When Meeting Heterophily[J]. arXiv preprint arXiv:2205.07308, 2022. ↩︎