图神经网络（GNN）入门笔记（1）——图信号处理与图傅里叶变换

一、信号处理：时域与频域

时域（Time Domain）是我们生活中常见的信号表示方式，以横轴为时间，纵轴为信号该时刻的强度（幅度），就可以得到一张时域图。打个比方，通过每时每刻的震动强度可以得到音频的时域图，或者画一张每小时的车流量，或者每日的降水量的图，也是时域图。
频域（Frequency Domain）则是一种数学上描述信号的方式。在频域中，唯一存在的波形是正弦波（Sine Wave，此处指代余弦和正弦这样的简谐波总称），通过不断叠加不同周期的正弦波，我们认为我们可以拟合所有的时域信号。画成图的话，频域是横轴为频率，纵轴为幅度的一个图形。

那么我们为什么需要频域呢？一方面，将复杂信号拆分为不同频率的正弦波，可以更简洁地对其进行表示。另一方面，通过频率分离，我们可以做到一些滤波操作。例如，一段音频里既有男生在说话，也有女生在说话，我们想把男女生的声音分开，而女生的声音频率要比男生的高，于是就可以根据频率来划分。保留男生的低频声音的滤波叫做低通滤波，而保留女生的高频声音的滤波叫做高通滤波。

二、傅里叶变换

2.1 傅里叶变换的本质

傅里叶变换的总公式如下：
$$F(\omega )={\int }_{-\infty }^{+\infty }f(t)e^{-i\omega t}dt$$
其中$\omega$表示频率，而$t$表示时间，这个操作就是把信号从时域转到频域上。
这个可怕的式子到底是什么意思呢？在这里非常推荐3Blue1Brown的图解。
让我们回到高中学过，或者高中的你正在学的复数概念。复数是二维的数，在横轴的实数轴上，又添加了一条纵轴的虚数轴。而单位虚数$i=\sqrt{-1}$。
著名欧拉公式就是$e^{ix}=\cos x+i\sin x$，这个公式有一个重要的几何意义，那就是$e^{ix}$代表着在复平面上，从(1,0)出发，沿着圆心为原点，半径为1的单位圆逆时针走$x$个单位的圆周，所落到的点代表的复数。
由此可知：

• $h(t)=e^{it}$ 代表一个每秒走1个单位的圆周的，在圆周上逆时针运动的点的轨迹函数。
• $h(t)=e^{-it}$是顺时针运动的轨迹函数。
• $h(t)=e^{-i\omega t}$是每$\frac{1}{\omega }$秒走1单位圆周（或者说每秒走$\omega$单位圆周），在圆周上逆时针运动的点的轨迹函数。
• $h(t)=f(t)e^{-i\omega t}$，则相当于我们将时域图上，横轴为时间，纵轴为幅度的$f(t)$，画到了圆上。时间从轴变成了角度，幅度从轴变成了幅长（到原点的距离）。圆上的一圈，我们能画原时域图$\frac{1}{\omega }$秒的内容。
  
  对于一个正弦波，大部分情况下只要你无限地画下去，它在圆周上都会分布得比较均匀，除了……当$\frac{1}{\omega }$秒恰好是原正弦波的一个周期时。此时，时域图上波峰的位置会在圆的一侧，而波谷的位置会在另一侧，就像这样：
  
  那么这时，我们将这个图形上所有的点都加起来，也就是，找它们二维上的几何平均点（质心），就可以近似得到原函数在$\omega$这个频率下的正弦波有多大的强度。而这种正弦波在原函数中持续时间越长（画了越多圈上图的图案），得到的$F(\omega )={\int }_{-\infty }^{+\infty }f(t)e^{-i\omega t}dt$也会越大。由此，我们就得到了原时域信号在不同频率下的强度。
  同样，也存在频域到时域上的傅里叶逆变换：
  $$f(t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F(\omega )e^{i\omega t}d\omega$$

2.2 从线性代数角度看离散傅里叶变换

能够被存储下来的信号，往往没有真正的连续形式，故实际中用得更多的，是离散傅里叶变换。既然离散了，我们显然也没办法把所有频率都做一遍了，此时可以近似地用 2 k π n ( k ∈ { 0 , . . . , n − 1 } \frac{2k\pi}{n}(k \in \{0,...,n-1\} n2kπ​(k∈{0,...,n−1}这一套频率来近似地拟合，下面我们记记$W_n=e^{\frac{2\pi }{n}i}$。
假设我们在时域上有一组采样点$f(t_0),...,f(t_{n-1})$，离散傅里叶变换的公式为：
$$F(W_n^k)=\sum _{j=0}^{n-1}f(t_j)W_n^{-kj}$$
相应地，逆傅里叶变换的公式为：
$$f(t_j)=\frac{1}{n}\sum _{k=0}^{n-1}F(W_n^k)W_n^{kj}$$
如果把傅里叶变换写成矩阵形式，就是这样：
$$\left(\begin{array}{c}F(0)\\ ⋮\\ F(n-1)\end{array}\right)=\left(\begin{array}{ccc}{W}_n^{-0}& \cdots & W_n^{-(n-1)}\\ ⋮& \ddots & ⋮\\ W_n^{-(n-1)}& \cdots & W_n^{-(n-1{)}^2}\end{array}\right)\left(\begin{array}{c}f(t_0)\\ ⋮\\ f(t_{n-1})\end{array}\right)$$

中间那个巨大的矩阵，我们记作$U$，于是上式可以简写为$F=Uf$。
同样，我们可以把逆傅里叶变换写成$f=\frac{1}{n}VF$，其中$V_{ij}=W_n^{ij}$
考虑$(VU{)}_{ij}={\sum }_{k=0}^{n-1}{W}_n^{ik}{W}_n^{-kj}={\sum }_{k=0}^{n-1}{W}_n^{k(i-j)}$
当$i=j$时，显然和为$(VU{)}_{ii}=n$。而$i\ne j$时，根据等比数列求和可得$(VU{)}_{ij}=\frac{W_n^{n(i-j)}-1}{W_n^{(i-j)}-1}$，由$W_n^n=1$可知此时$(VU{)}_{ij}=0$。

三、拉普拉斯算子

拉普拉斯算子（Laplace Operator）写作一个正着的三角，表示梯度的散度，简单来说，对于$n$维函数（注意这里是$n$维而不是$n$个采样点）$f(x_1,...,x_n)$有：
$$\Delta f=\nabla \cdot \nabla f=\left(\begin{array}{ccc}\frac{\partial }{\partial x_1}& \dots & \frac{\partial }{\partial x_n}\end{array}\right)\left(\begin{array}{c}\frac{\partial }{\partial x_1}\\ ⋮\\ \frac{\partial }{\partial x_n}\end{array}\right)f=\sum _{i=1}^n\frac{{\partial }^{2}f}{\partial x_i^2}$$
那么拉普拉斯算子和频率有什么关系呢？写成连续形式或许不够直观，我们可以考虑一个二维离散的场景，即，对一张灰度图使用拉普拉斯算子。函数$f(x,y)$表示的是灰度图上$(x,y)$点的灰度。在此情景下，我们近似地将导数$\frac{\partial f}{\partial x}={\lim }_{h\to 0}(f(x+h)-f(x))$写作$f(x+1)-f(x)$或$f(x)-f(x-1)$，则：
$\frac{{\partial }^{2}f(x_0,y_0)}{\partial x_0^2}+\frac{{\partial }^{2}f(x_0,y_0)}{\partial y_0^2}$
$=\frac{\partial f(x_0+1,y_0)}{\partial x_0}-\frac{\partial f(x_0,y_0)}{\partial x_0}+\frac{\partial f(x_0,y_0+1)}{\partial y_0}-\frac{\partial f(x_0,y_0)}{\partial y_0}$
$=f(x_0+1,y_0)+f(x_0,y_0+1)+f(x_0-1,y_0)+f(x_0,y_0-1)-4f(x_0,y_0)$
而这反映了这个像素点处灰度值的变化大小，某种意义上相当于此处的频率。这也是CV中的拉普拉斯卷积核的来历：
$$\left(\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right)$$
既然拉普拉斯算子是干这个用的，可以想见它和傅里叶变换会有些关系。为了寻求这个关系，我们可以回顾一下线性代数中特征值和特征向量的概念：

特征值和特征向量
概念：对于一个矩阵$A$，若存在向量$u_i$和常数值${\lambda }_i$，使得$A{u}_i={\lambda }_i{u}_i$，则称${\lambda }_i$为$A$的特征值，$u_i$为$A$的特征向量。
求法：则$(A-{\lambda }_{i}I)u_i=0$，通过解$|A-\lambda I|=0$可以得到$n$组解，再相应带入回原式可以得到相应的${\lambda }_i$。矩阵的特征值集合是唯一的，而每个特征值对应的特征向量构成一个特征子空间。
矩阵的对角化：若$A$有$n$个线性无关的特征向量$U=(u_1,...,u_n)$，则$A=Udiag({\lambda }_1,...,{\lambda }_n)U^{-1}$
利用对角化计算矩阵向量乘：将向量$x$映射到空间$U$中，得到$x=x_1{u}_1+...+x_n{u}_n$，则$Ax=x_1{\lambda }_1{u}_1+...+x_n{\lambda }_n{u}_n$

那么，假如有这样一个函数$f$和常数值$\lambda$，使得$\Delta f=\lambda f$，岂不是这就是$\Delta$的特征值和特征函数了吗？那么再想想，什么样的函数，经过两次求导，还能等于一个常数乘以它本身呢？答案是……$\sin x$，$\cos x$和$e^x$，对吧？事实上，$\frac{d^2{e}^{ix}}{d{x}^2}{e}^{iax}=-a^2{e}^{iax}$，也即$e^{iax}$系列函数是拉普拉斯算子的特征函数。
我们可不可以像利用特征向量计算矩阵向量乘一样，利用这特征函数来计算拉普拉斯算子呢？那么第一步，我们就需要把一个任意函数$f(x)$，映射到$e^{iax}$的函数族空间中。实际上，套一下傅里叶函数的公式可知：
$$f(x)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }({\int }_{-\infty }^{+\infty }f(t)e^{-i\omega t}dt)e^{i\omega x}d\omega$$
记$c(\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }f(t)e^{-i\omega t}dt$，如果它收敛到一个常数值，则这个常数值是和$f$整体有关，而和自变量$x$的具体取值无关的。$f(x)$也就就此被映射到了特征函数族中。这时候，想做拉普拉斯算子，就是：

$$\Delta f(x)={\int }_{-\infty }^{+\infty }{\omega }^{2}c(\omega )e^{i\omega x}d\omega$$

四、图的拉普拉斯矩阵

设$A$为图的邻接矩阵（即当i和j存在边时$A_{ij}=1$，否则$A_{ij}=0$，在此我们只考虑无向图，即$A$为对称矩阵），而$D$为图的度数矩阵（对角阵，$D_{ii}$表示点$i$有多少条邻边），定义图拉普拉斯矩阵为：
$$L=D-A$$
形象化地理解，我们考虑这样一个实际问题。有一个公园建在山上，公园内有$n$个景点，一些景点之间有路相连。我们想要分析这个公园某个景点有多陡峭。应用拉普拉斯矩阵的话，是这样的：
$$\Delta h(u)=\sum _{v\in Neighbor(u)}(h(u)-h(v))$$
$$\Delta h=\left(\begin{array}{c}\Delta h(1)\\ ⋮\\ \Delta h(n)\end{array}\right)=(D-A)h=LA$$
当然，就这么直接用似乎还有些问题，直接用上面这个$L$的话，似乎一个节点邻居越多，最后得到的变化值就越大，因此实际中我们一般会对其进行归一化，比如这样：
$$\hat{L}=D^{-1}L=I-D^{-1}A$$
不过这种非对称归一化会有点小问题，我们假设点1有很多邻居，而点2只有点1一个邻居。在归一化时，$\Delta h(1)$只除了点1的度数，而$\Delta h(2)$只除了点2的度数。在一些问题（例如图上标签预测问题）中，我们认为此时2受到1的影响过高了。解决方案则是，在归一化过程中，点2受到点1的影响，不仅要除以和点2度数相关的量，也要除以和点1度数相关的量。具体来说，就是这样做对称归一化：
$$\hat{L}=D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}=I-D^{-\frac{1}{2}}A{D}^{-\frac{1}{2}}$$

五、图傅里叶变换

一维信号的傅里叶变换，本质是用不同频率的正弦波来拟合原信号。那我们应该如何构造一组图上的基本信号，来拟合原图的函数呢。
考虑把信号记为$f=(f(1),...,f(n){)}^T$，然后考虑用下式，通过引入平方消除正负性的影响，来反映整张图的波动性（variation）：
$$f^{t}Lf=\sum _{i,j,i<j}{A}_{ij}(f(i)-f(j){)}^2$$
首先我们可以发现，对于$u_1=(\frac{1}{\sqrt{n}},...,\frac{1}{\sqrt{n}}{)}^T$（这里对向量的模做了归一化），$L{u}_1=0$，即$u_1$是$L$的一个特征向量，$0$为其特征值，而且它能够使得原图的波动最小。接下来，我们不断重复以下过程：
$$u_k=\underset{||f||=1,f|u_1,..u_{k-1}}{\min }{f}^{t}Lf$$
我们就会得到一组$L$的特征向量正交基，它们也是一组使得原图波动最小的基本图信号，而且对应特征值${\lambda }_k=u_k^{t}L{u}_k$。
设$\tilde{x}$为图信号$x$在频域下的表示，则$x=u_1{\tilde{x}}_1+...+u_n{\tilde{x}}_n=U\tilde{x}$。相应的，我们可以定义将$x$表示在频域上的图傅里叶变换为$\tilde{x}=U^{T}x$，而特征值$({\lambda }_1,...{\lambda }_n)$，就是$n$组频率。
下图就是一个不同频率上拉普拉斯特征向量的可视化表现：

六、图滤波与图卷积

我们之前提过，在声音信号中，通过改变原信号不同频率的强弱，可以做到诸如分离男女声音之类的功能，这个过程就叫滤波。此外还有一个经常提到的概念叫卷积，假设$f$为原函数，$g$为卷积核，$f$和$g$之间做卷积写成公式是这样的：
$$(f\ast g)(t)={\int }_{-\infty }^{\infty }f(t^{\prime })g(t-t^{\prime })d{t}^{\prime }$$
而我们有：

卷积定理：若$\hat{f}$是$f$在频域上的表示，则$\widehat{(g{\ast }_{g}f)}(t)=(\hat{f}\cdot \hat{g})(\omega )={\int }_{-\infty }^{+\infty }\hat{f}(\omega )\hat{g}(\omega )d\omega$。

虽然在这里我不打算展开证明，但可以简单介绍一下，利用这个定理，我们可以通过（1）将$f$和$g$用$O(n\log n)$从时域转化为频域 （2）对二者在频域进行$O(n)$点积 （3）频域$O(n\log n)$转化回时域 的三个步骤，将本来需要$O(n^2)$计算的卷积在$O(n\log n)$中快速完成。
让我们放眼第二步，你发现了什么？——没错，在频域上进行点积，不就是对$f$在不同频率的强度进行调整吗？所以卷积和滤波本质上是相同的操作。
好，因为我们已经定义过图傅里叶变换了，所以可以以同样流程考虑图的情况，设$f$为图结点的特征组成的向量，$g$是卷积核，则：
$$(g{\ast }_{g}f)=U(U^{T}g\cdot U^{T}f)$$
令$W=diag(U^{T}g)$，则$(g{\ast }_{g}f)=UW{U}^{T}f$
也即，如果我们希望对图信号进行处理，只需要找到（使用各种学习方法）合适的一组$W$的$n$个参数即可。这就是谱域（Spectral Domain）上的图卷积的基本思想。
当然，到目前为止，我们得到的这个信号处理方式想要实际应用还存在一些问题：

1. $W$是和图形状强相关的，参数量和图的点数也有关，难以适用于不同的图。
2. 需要计算图傅里叶算子（拉普拉斯矩阵的特征向量），这个过程是$O(n^3)$的，$n$为图的点数，复杂度较高且和图规模相关。
3. 只能用于无向图（需要拉普拉斯矩阵对称以保证可以找到正交基特征向量）。
4. 难以对局部空间信息进行处理（毕竟所有频率都是全局的）。
   我们将在下一节继续讨论这些问题。