图拉普拉斯的前世今生

2021-06-04

字数统计: 5.7k | 阅读时长≈ 23 分钟

图拉普拉斯的前世今生博客搬运工 # 1. 拉普拉斯算子

1.1 散度

二维平面上的通量

先看一道物理题：小明乘帆船出行，刮来一阵妖风，假设帆的面积为 $S$ , 和妖风的夹角为 $\alpha$ ，妖风在每单位面积上的垂直风压为 $S$ ，求妖风对帆的推动力。

显然答案是 $F=S\times P\times sin(\alpha)$ 。如果用 $\overrightarrow{A}$ 表示风的向量 $||\overrightarrow{A}||=P$ ，用 $\overrightarrow n$ 表示帆的法向量，那么 $\overrightarrow n$ 与风的夹角即为 $\frac{\pi}{2}-\alpha$ ，那么此时 $F=S·\overrightarrow A·\overrightarrow n$ 。

三维曲面上的通量

假设题目中的帆是空间中的一个曲面 $\sum$ ，在 $\sum$ 中每一点的法向量用 $\overrightarrow n$ 表示，空间中的风力用 $\overrightarrow A$ 表示，那么此时 $F=\displaystyle\iint_\sum\overrightarrow A·\overrightarrow ndS$ ，这就是通量 $\Phi$ 的定义 $\Phi_A(\sum)=\displaystyle\iint_\sum\overrightarrow A·\overrightarrow ndS$ 。
封闭曲面的通量

上图中红色箭头表示曲面在此处的法向量，在计算通量的时候我们通常认为向量场会穿过曲面，不会被挡住，此时 $\Phi_A(\sum)= \displaystyle∮_\sum\overrightarrow A·\overrightarrow ndS$ 。容易看出，根据向量乘法原理，上图中射入曲面和射出曲面的通量一正一负，相互抵消了，最终整个封闭曲面的总通量为0。（其实通量的概念和高中所学的磁通量一样，磁场中，通过封闭曲面的磁通量为0）

散度

所谓散度，其实从名字上来看描述的是：向量场中某个点它的发散程度，向量在这一点是向外发散的，还是向内收敛的，它的发散或收敛强度如何。如下图所示：

显然，这一向量场中，红色曲面的总统量为负，而绿色曲面上的总通量为正。如果我们无限缩小这个曲面，直至接近一个点 $\{x\}$ ，用 $\delta V$ 表示这个曲面围成的体积，那么：
$\Delta A=\mathop{lim}\limits_{\delta V\rightarrow x}\frac{\Phi_A(\sum)}{\delta V}$
$\Delta A$ 就表示点 $\{x\}$ 处的散度。上图中，红圈的圆心处的散度为负，而绿圈的圆心出散度为正。

因此，散度衡量了一个点处的向量场是被发射还是被吸收，或者说，对于散度为正的点，散度越大，意味着相应的向量场越强烈地在此发散，而对于散度为负的点，意味着相应的向量场在此汇聚。

1.2 拉普拉斯算子

定义

拉普拉斯算子是n维欧式空间中的一个二阶微分算子，定义为梯度 $\Delta f$ 的散度 $\Delta·\Delta f$ 。
涵义

以三维为例，在直角坐标系中，一个函数 $f(x,y,z)$ 在点 $(x_0,y_0,z_0)$ 处的梯度是一个向量 $(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z})|_{x=x_0,y=y_0,z=z_0}$ 。因此函数f的梯度函数为 $\Delta f=\frac{\partial f}{\partial x}·\overrightarrow i+\frac{\partial f}{\partial y}·\overrightarrow j+\frac{\partial f}{\partial z}·\overrightarrow k$ ，这就构成了一个在三维空间下的向量场。我们再对这个向量场 $\Delta f$ 求散度 $\Delta·\Delta f$ 即可得到 $f$ 的拉普拉斯算子 $\Delta^2 f$ 。
这样定义的原因

让我们想像一座山，根据梯度的定义，在山峰周围，所有的梯度向量向此汇聚，所以每个山峰处的拉普拉斯算子为负；而在山谷周围，所有梯度从此发散，所以每个山谷处的拉普拉斯算子为正。所以说，对于一个函数，拉普拉斯算子实际上衡量了在空间中的每一点处，该函数梯度是倾向于增加还是减少 。

（调和函数， $\Delta^2f=0$ ）

2. 图拉普拉斯

2.1 图论下的函数

函数即从定义域到值域的一个映射。在图中，如果用 $V$ 表示图中所有顶点构成的结合， $E$ 表示所有边构成的集合，我们可以定义一个图函数 $F_G:V\rightarrow R$ ，使得图中每个节点 $v\in V$ ，都被映射到 $R$ 中的一个实数。以下图社交网络为例：

图中每一条边表示两个人之间信息的流通程度。现在我们想要分析这个社交网络上的信息传播，我们不仅需要知道信息流通的程度，我们还要知道每个人发动态的活跃程度，于是我们现在给这个图一个函数 $F_G$ ，使得：

F_G(A)=3,\ F_G(B)=5,\ F_G(C)=-2,\ F_G(D)=1,\ F_G(E)=-5

这里的负数可以理解为：C和E是聊天终结者，可以阻止信息的传播。这样我们可以得到下面这张图：

2.2 图函数的梯度

梯度描述的时函数在每一点处正交方向上的变化（笛卡尔坐标系中的坐标轴），如 $f(x,y,z)$ 的梯度在x方向的分量为 $\frac{\partial f}{\partial x}=\frac{f(x+dx)-f(x)}{(x+dx)-x}$ 。类比到图中，我们把节点看做欧式空间中的一个点，把边看做是正交的两个分量（笛卡尔坐标系中的坐标轴），并且对于两个节点，我们定义其距离d是其边权值的倒数。那么对于一个节点 $v_0$ ，我们认为其梯度在一条通向 $v_1$ 的边 $e_{01}$ 上的分量为：

\frac{(F_G(v_0)-F_G(v_1))}{d_{01}}=(F_G(v_0)-F_G(v_1))·e_{01}

为了计算梯度，我们首先定义下面的矩阵 $K_G$ ：

其中每一行代表一个点，每一列代表一条边，使得对于每条边，如果该条边冲改点发出的，且权值为X，那么将矩阵中对应的位置的元素置为 $-X$ ，如果该条边指向改点，那么将对应元素设为 $+X$ 。

同时根据之前图函数 $F(G)$ 的定义，我们有下面的列向量 $f_G$ ：

下面我们计算 $K_G^T\times f_G$ ：

经观察可以发现，计算结果即为根据 $(F_G(v_0)-F_G(v_1))·e_{01}$ 计算出的整个图的梯度向量，其中每一行表示梯度在这一条边上的分量。因此，对于图G，我们有 $\Delta_G=K_G^T$ ，使得 $\Delta_GF_G=K_G^T\times f_G$ 。

2.3 拉普拉斯矩阵

2.3.1 拉普拉斯算子

在函数中，拉普拉斯算子定义为函数梯度的散度，即每一点上其梯度的增加/减少。类别到图中，我们显然可以这样定义图函数的散度：即从该节点射出的梯度，减去射入该节点的梯度。我们可以按照下面方法计算散度，即将原来的梯度向量左乘一个矩阵 $K_G'$ ：

$K_G'$ 中每一行代表一个点，每一列代表一条边，使得对于每个点每条边，如果该条边从该点发射出去，则将矩阵中对应的这一元素置为-1，如果该条边指向该点，则将对应的元素置为 +1。将 $K_G$ 中正元素变成1，负元素变成-1即可得到 $K_G'$ 。

这样我们就有：

\begin{align} \Delta_G^2F_G&=K_G'\times\Delta_GF_G\\ &=K_G'\times(K_G^T\times f_G)\\ &=(K_G'K_G^T)f_G \end{align}

于是我们就得到了图函数中拉普拉斯算子的定义 $\Delta^2=K_G'K_G^T$ ，即拉普拉斯矩阵。需要注意的是，拉普拉斯矩阵的值与图中每一条边的方向无关，所以拉普拉斯矩阵一般用于表述无向图。计算 $K_G'K_G^T$ 的值，我们可以得到：

计算结果为对称矩阵，对角线为每个节点的度，而其余元素则是负的邻接矩阵，于是可以得到拉普拉斯的经典算式：

L=D-W

拉普拉斯矩阵的一个重要性质就是：具有n个非负特征值 $\lambda_1,\lambda_2,...\lambda_n$ ，且 $0=\lambda_1\leq\lambda_2\leq...\leq\lambda_n$ 。同时根据拉普拉斯矩阵和瑞利熵的概念 $R(L,h)=\frac{h^*Lh}{h^*h}$ ，通过拉格朗日乘子法可以得出：瑞利熵的最大值，等于L的最大特征值，瑞利熵的最小值等于L的最小值。

2.3.2 差分角度看拉普拉斯

在欧式空间的函数中，拉普拉斯算子是这样定义的： $\Delta=\sum\limits_i\frac{\delta^2}{\delta x_i^2}$ ，即非混合二阶偏导的和，其余所有的XXXX拉普拉斯算子都是其中的一个特例，或者是其某种情况下的一种近似。

图像上的拉普拉斯

图像是一种离散数据，那么我们图像上的拉普拉斯必然要进行离散化处理：
- 导数的定义： $f'(x)=\frac{\delta f(x)}{\delta x}\approx f(x+1)-f(x)$
- 二阶导数： $\frac{\delta^2f(x)}{\delta x^2}=f''(x)\approx f'(x)-f'(x-1)=f(x+1)+f(x-1)-2f(x)$
这样我们可以得到两个结论：
- 结论1：二阶导数近似等于其二阶差分
- 结论2：二阶导数等于其在所有自由度上微扰后获得的增益。（自由度可以理解为笛卡尔坐标系中的坐标轴，比如1维函数的自由度可以理解为2，有+1和-1两个方向，而2为函数的自由度则为4，分别是x+1，x-1，y+1，y-1四个方向）
如上，对于二维图像来说共有两个方向，4个自由度可以变化，即如果对 $(x,y)$ 处的像素进行扰动，其可以变成四种状态 $(x+1,y),(x-1,y),(x,y+1),(x,y-1)$ 。当然如果我们考虑对角线方向也是一种自由度的话，还会增加四种，这个还是看个人的设定。如果我们只考虑二维图像上4个自由度的变化，那么可以得到：

\begin{align} \Delta&=\frac{\delta^2f(x,y)}{\delta x^2}+\frac{\delta^2f(x,y)}{\delta y^2}\\ &\approx f(x+1,y)+f(x-1,y)-2f(x,y)+[f(x,y+1)+f(x,y-1)-2f(x,y)]\\ &= f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y) \end{align}

上式可以理解为对 $(x,y)$ 点的像素进行轻微扰动后，使其变化到相邻像素后得到的增益。

这给我们一种形象的结论：拉普拉斯算子就是在所有自由度上进行微小变化后获得的增益。
图上的拉普拉斯

推广到图上，对于N个节点的图，假设其邻接矩阵为W，这个图的自由度是多少呢？最多是N。如果这个图是完全图，那么对任意一个节点进行微扰，它就可能变成任意一个节点。此时类比图像中的 $f(x,y)$ ，图中的函数应该就是一个N维向量 $f=(f_1,...,f_N)$ ，其中 $f_i$ 表示函数在节点 $i$ 出的函数值。

根据图像上的拉普拉斯，拉普拉斯算子就是在所有自由度上进行微小变化后获得的增益。对于Graph，从节点i变化到节点j增益是多少呢？比较容易想到的就是利用节点之间边的权值来计算：
$\sum\limits_{j\in\mathcal N_i}W_{ij}[f_j-f_i]$
因此，对于图来说，其拉普拉斯算计如下：
$(\Delta f)_i=\sum\limits_i\frac{\delta^2f}{\delta i^2}\approx\sum\limits_{j\in\mathcal N_i}W_{ij}[f_j-f_i]$
上式中 $j\in\mathcal N_i$ 可以去掉，因为节点i和j不相连的话， $W_{ij}=0$ ，所以上式可以继续化简为：
$\begin{align} \sum\limits_{j\in\mathcal N_i}W_{ij}[f_j-f_i]&=\sum\limits_jW_{ij}f_j-\sum\limits_jW_{ij}f_i\\ &=(Wf)_i-(Df)i\\ &=[(W-D)f]_i \end{align}$
上式对于任何i都成立，也就有： $\Delta f=(W-D)f$ ，即图上的拉普拉斯算子为 $W-D$ 。

2.3.3 瑞利熵

普通瑞利熵

$R(A,x)=\frac{x^TAx}{x^Tx}$

其中A为n阶Hermitian阵，在机器学习中通常处理的都是实数，所以A为实对称矩阵（图中的拉普拉斯是不是满足A的条件？）。通过拉格朗日乘子可以证明：
- R的最大值就是A的最大特征值，R的最小值就是A的最小特征值
- x的解，就是A的特征值对应的特征向量
瑞利熵经常出现在聚类和降维任务中，因为降维、聚类任务往往会导出最大化最小化瑞利熵的式子，从而通过特征值分解的方式找到降维空间。

PCA、谱聚类、RatioCut等算法都是瑞利熵的一个应用。
泛化瑞利熵

对于 $R(L,x)=\frac{x^TLx}{x^Tx}$ ，其分子是一般形式，而分母是 $x^TDx$ 在D取单位阵时的特殊情况，将其一般化：
$R(L,h)=\frac{x^TLx}{x^TDx}$
同理，通过拉格朗日乘子可以导出：
$R(L,h)=\frac{x^T[D^{-\frac{1}{2}}LD^{-\frac{1}{2}}]x}{x^Tx}$
这样又回到了普通瑞利熵问题了。（ $D^{-\frac{1}{2}}LD^{-\frac{1}{2}}$ 这个和正则化的图拉普拉斯是不是一样？）

Fisher线性判别器（LDA）、Ncut等算法就是泛化瑞利熵的一种应用。

3. 谱聚类

（该部分主要从刘建平博文https://www.cnblogs.com/pinard/p/6221564.html搬运而来）

谱聚类是从图论中演化出来的算法，后来在聚类中的到了广泛的应用。它的主要思想是：把所有数据看作空间中的点，这些点可以用边连接起来，距离较远的点之间的边的权重较低，距离近的边的权重较高，这样就构成了图。通过对所有数据点组成的图进行切图，让切图后得到的子图之间的边的权重之和尽可能低，而子图内边权重之和尽可能高，从而达到聚类的目的。

3.1 构建无向权重图

对于图G，我们通常用顶点集合V和边集合E来表述，即 $G(V,E)$ 。对于V中的任意两个顶点，可以有边连接，也可以没有边连接。我们定义 $w_{ij}$ 表示点 $v_i$ 和 $v_j$ 之间的权重。因为我们的图是无向的，所以 $w_{ij}=w_{ji}$ 。

对于有边连接的两个点 $v_i$ 和 $v_j$ ， $w_{ij}>0$ 。对于没有边连接的两个点 $v_i$ 和 $v_j$ ， $w_{ij}=0$ 。对于图中的任意一个点 $v_i$ ，它的度 $d_i$ 定义为和它相连的所有边的权重之和，即：

d_i=\sum\limits_{j=1}^nw_{ij}

这样我们可以得到一个n阶对角矩阵D：

批注 2020-12-06 142153.png

它只有对角线元素有值，对应第i行的第i个点的度。利用所有点之间的权重值，我们可以得到图的邻接矩阵W，它是一个n阶对称矩阵，第i行第j个元素表示点 $v_i$ 和点 $v_j$ 之间边的权重。

另外，对于点集V的一个子集 $A\subset V$ ，我们定义：

|A|:=子集A中节点的个数\\ vol(A):=\sum\limits_{i\in A}d_i

3.2 构造相似矩阵

在构造加权无向图的过程中，需要使用到边的权重来构造邻接矩阵。但是我们只有数据点的定义，并没有边的定义，如何构造图的邻接矩阵呢？

一个基本的思想就是：距离较远的两个点之间权重值较低，而距离较近的两个点之间的边的权重较高。通常我们使用样本点距离度量的相似矩阵S来获得邻接矩阵W。常用的构建邻接矩阵W的方法有三类： $\epsilon$ -邻近法，K邻近法和全连接法。

$\epsilon$ -邻近法

设置一个距离阈值 $\epsilon$ ，然后用欧式距离 $s_{ij}=||x_i-x_j||_2^2$ 度量任意两个点 $x_i$ 和 $x_j$ 的距离，然后根据 $s_{ij}$ 和 $\epsilon$ 的大小关系来定义邻接矩阵W如下：
$\begin{equation} w_{ij}= \left\{ \begin{array}{lr} 0&\quad s_{ij}>\epsilon\\ \epsilon&\quad s_{ij}\leq\epsilon \end{array} \right. \end{equation}$
两点之间的权重要么是 $\epsilon$ ，要么是0，没有其他信息了。由于距离远近度量很不精确，因此在实际应用中，我们很少使用 $\epsilon$ -近邻法。
K近邻

使用KNN算法遍历所有的样本点，取每个样本的最近的k个点作为近邻，只有和距离最近的k个点之间的 $w_{ij}>0$ ，其余点之间的 $w_{ij}=0$ 。但是这种方法会导致得到的邻接矩阵W非对称，为了解决这种问题一般采取下面两种方法之一：
- 第一种K近邻法只要一个点在另一个点的K近邻中，就保留 $s_{ij}$
  $\begin{equation} w_{ij}=w_{ji} \left\{ \begin{array}{lr} 0& x_i\notin KNN(x_j)\ and\ x_j\notin KNN(x_i)\\ exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& x_i\in KNN(x_j)\ or\ x_j\in KNN(x_i) \end{array} \right. \end{equation}$
- 第二种K邻近法是必须两个点互为K近邻，才能保留 $s_{ij}$
  $\begin{equation} w_{ij}=w_{ji} \left\{ \begin{array}{lr} 0& x_i\notin KNN(x_j)\ or\ x_j\notin KNN(x_i)\\ exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& x_i\in KNN(x_j)\ and\ x_j\in KNN(x_i) \end{array} \right. \end{equation}$
全连接法

和前两种方法相比，这种方法所有点之间的权重值都大于0，因此称之为全连接法。我们可以选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和Sigmoid核函数。最常用的时搞死核函数RBF，此时相似矩阵和邻接矩阵相同：
$w_{ij}=s_{ij}=exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})$

在实际应用中，使用全连接方法来建立邻接矩阵是最普遍的，而在全脸接法中使用高斯径向核RBF是最普遍的。

3.3 无向图切图

对于无向图G的切图，我们的目标是将图 $G(V,E)$ 切成相互没有连接的k个子图，每个子图中点的集合为： $A_1,A_2,...,A_k$ ，它们满足 $A_i\cap A_j=\O$ ，且 $A_1\cup A_2\cup...\cup A_k=V$ 。对于任意两个子图 $A,B\subset V,A\cap B=\O$ ，我们定义A和B之间的切图权重为：

W(A,B)=\sum\limits_{i\in A,j\in B}w_{ij}

那么我们定义切图cut为：

cut(A_1,A_2,...,A_k)=\frac{1}{2}\sum\limits_{i=1}^kW(A_i,\bar A_i)

其中 $\bar A_i$ 表示 $A_i$ 的补集，即除 $A_i$ 外的其他V的子集的并集。那么如何切图可以让子图内的点权重和高，子图间的权重和低呢？一个自然的想法就是最小化 $cut(A_1,A_2,...,A_k)$ ，但是发现，这种最小化方式存在下图所示的问题：

批注 2020-12-06 145842.png

此时，如果我们选择一个权重最小的边缘点，比如C和H之间进行cut，只要可以最小化 $cut(A_1,A_2,...,A_k)$ ，但是得到的却不是最优的切图。如何避免这种切图，并且找到类似图中“Best Cut”这样的最优切图呢？我们需要对每个子图的规模做出限定，下面就来介绍两种常用的切图方式。

3.3.1 RatioCut切图

RatioCut切图为了避免上图中的最小切图，对每个切图，不光考虑最小化 $cut(A_1,A_2,...,A_k)$ ，它还同时考虑最大化每个子图的数量，即：

RatioCut(A_1,A_2,...,A_k)=\frac{1}{2}\sum\limits_{i=1}^k\frac{W(A_i,\bar A_i)}{|A_i|}

那么如何最小化上面这个RatioCut函数呢？牛人们发现，RatioCut函数可以通过如下方式表示：

3.3.2 Ncut切图

4. GCN

5. 参考链接

打赏

版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！