图/超图拉普拉斯分解 Graph/Hypergraph Laplacian eigendecomposition

大家好，还是之前的风格，简单粗暴，直接开干！

一、图拉普拉斯特征分解

1. 数学公式

直接上公式，下边公式是图的归一化 拉普拉斯特征分解公式：
$$\Delta =I-D^{-\frac{1}{2}}A{D}^{-\frac{1}{2}}=U^T\Lambda U$$

各参数的含义是，给定一个图G，假如它有n个节点和m个边:

• $\Delta$ ：拉普拉斯特征矩阵，我们要先计算出这个矩阵，才方便进行特征分解， 矩阵shape是 $n\times n$

• $I$：单位矩阵， 矩阵shape是 $n\times n$

• $D$：节点度矩阵，矩阵shape是 $n\times n$

• $A$ ：邻接矩阵，矩阵shape是 $n\times n$

• $U$：特征向量， 矩阵shape是 $n\times n$，用代码获取的一般是按照对应特征值从小到大的顺序排列的，每一行对应每个节点的特征值，每一列对应一个特征值。对应特征值小的特征向量包含了更多图整体的结构信息，对应特征值大的特征向量包含了更多这个节点附近的局部信息。

• $\Lambda$：特征值， 矩阵shape是 $1\times n$，一般从函数获得后都是自动排好序的，从小到大

2. 参考代码

参考代码 Github

下边的代码是我加了一些个人注释的，便于理解：

import torch
import torch.utils.data
import numpy as np


def eig(sym_mat):
	# 特征分解函数
    # (sorted) eigenvectors with numpy
    EigVal, EigVec = np.linalg.eigh(sym_mat) # 特征分解函数，输入图邻接矩阵A，返回特征值和特征向量

    # for eigval, take abs because numpy sometimes computes the first eigenvalue approaching 0 from the negative
    eigvec = torch.from_numpy(EigVec).float()  # [N, N (channels)] # 转一下type，无伤大雅
    eigval = torch.from_numpy(np.sort(np.abs(np.real(EigVal)))).float()  # [N (channels),] #其实分解函数返回的时候就已经是排过序的了
    return eigvec, eigval  # [N, N (channels)]  [N (channels),]


def lap_eig(dense_adj, number_of_nodes, in_degree):
    """
    拉普拉斯特征矩阵函数
    Graph positional encoding v/ Laplacian eigenvectors
    https://github.com/DevinKreuzer/SAN/blob/main/data/molecules.py
    """
    dense_adj = dense_adj.detach().float().numpy() # 邻接矩阵
    
    # 节点入度，这里in_degree输入的是每个节点的入度，现在还是一个一维的向量
    # 大概长这样：[2, 1, 4, 2, ..., 3, 3]，长度就是number_of_nodes，对于无向图来说也是出度矩阵
    in_degree = in_degree.detach().float().numpy() 

    # Laplacian
    A = dense_adj # 邻接矩阵
    N = np.diag(in_degree.clip(1) ** -0.5) # diag就是把 1xn 的向量转为 nxn 的矩阵
    L = np.eye(number_of_nodes) - N @ A @ N #这就是上文说的计算拉普拉斯特征矩阵的公式

    eigvec, eigval = eig(L) # 上边那个特征分解函数
    return eigvec, eigval  # [N, N (channels)]  [N (channels),]

二、超图拉普拉斯特征分解

1. 数学公式

下边公式是超图的归一化 拉普拉斯特征分解公式:

$$L=I-D_v^{-\frac{1}{2}}HW{D}_e^{-1}{H}^T{D}_v^{-\frac{1}{2}}=U^T\Lambda U$$

各参数的含义是，给定一个超图G，假如它有n个节点和m个边:

• $L$ ：拉超图普拉斯特征矩阵，我们要先计算出这个矩阵，才方便进行特征分解，矩阵shape是 $n\times n$

• $I$：单位矩阵， 矩阵shape是 $n\times n$

• $D_v$：节点度矩阵，矩阵shape是 $n\times n$

• $D_e$：超边度矩阵，矩阵shape是 $m\times m$

• $H$ ：关联矩阵，矩阵shape是 $n\times m$

• $W$：包含超边权重的对角矩阵，如果没有明确给出超边权重矩阵，一般超边的权重矩阵就按超边的度矩阵算，矩阵shape是 $m\times m$

• $U$：特征向量， 矩阵shape是 $n\times n$

• $\Lambda$：特征值， 矩阵shape是 $1\times n$

2. 参考代码

参考代码 Github

这是我自己写的参考代码，由于接收的数据不太理想，所以代码中可能含有除公式外一些额外的东西：

import numpy as np

def getEigvec(HT):
    H, Dv, De, W = getHypergraph(HT)
    eigval, eigvec = eig(H, Dv, De, W)

    return eigval, eigvec

def getHypergraph(HT):
    """
    This function aims to convert a graph to a hypergraph.
    
    graph: 
        HT: a hypergraph incidence matrix, shape is [num_hyperedges, num_hypernodes], contains 0 and 1.
    
    return Hypergraph:
        H: a hypergraph incidence matrix, shape is [num_hypernodes, num_hyperedges], contains 0 and 1.
        Dv: node degree matrix --> [num_nodes, num_nodes]
        De: hyperedge degree matrix --> [hyperedge_num, hyperedge_num]
        W: hyperedge weight matrix --> [hyperedge_num, hyperedge_num]
        
    """
    H = HT.T

    De = np.diag(np.sum(H, axis = 0)) # hyperedge degree matrix --> [num_hyperedges, num_hyperedge]
    Dv = np.diag(np.sum(H, axis = 1)) # nodes degree matrix  --> [num_nodes, num_nodes]
    W = De    # hyperedge weight matrix --> [hyperedge_num, hyperedge_num]
    
    return H, Dv, De, W

def eig(H, Dv, De, W):
    """
    This function aims to get hypergraph eigendecomposition.
    Reference: 
        https://github.com/DevinKreuzer/SAN/blob/main/data/molecules.py
        https://github.com/jw9730/tokengt/blob/main/large-scale-regression/tokengt/data/algos.py
    
    H: a hypergraph incidence matrix, shape is [num_hypernodes, num_hyperedges], contains 0 and 1.
    Dv: node degree matrix --> [num_nodes, num_nodes]
    De: hyperedge degree matrix --> [hyperedge_num, hyperedge_num]
    W: hyperedge weight matrix --> [hyperedge_num, hyperedge_num]
    
    return:
        EigVal: eigenvalue --> [num_nodes, 1]
        EigVec: eigenvector --> [num_nodes, num_nodes]
    """
    
    Dv_half = np.diag(Dv.sum(axis = 0).clip(1) ** -0.5)
    De_1 = np.diag(De.sum(axis = 1).clip(1) ** -1)
    num_nodes = len(Dv)
    
    delta = np.eye(num_nodes) - Dv_half @ H @  De_1 @ H.T @ Dv_half
    
    EigVal, EigVec = np.linalg.eigh(delta)
    # EigVal = np.sort(np.abs(np.real(EigVal)))
    idx = EigVal.argsort()
    EigVal, EigVec = EigVal[idx], np.real(EigVec[:, idx])


    return EigVal, EigVec