Heat Kernel Based Community Detection

heat kernel是一种局部扩散方法，所谓的局部扩散方法就是说给定一个种子节点，我以这个种子节点为源，向周围扩散，就可以找到这个种子节点所属于的社区。相似的思想有pagerank（Local Graph Partitioning using PageRank Vectors）。heat kernel方法相比于pagerank在于它将更大的权重放在了距离种子节点较近的节点，收敛的更快，在真实网络中证实了这种性质是非常有用的。

PRELIMINARIES

我们给定图$G=(V,E)$，有n个节点，$A$是邻接矩阵，因为G是无向图，所以A是对称矩阵。令D是度矩阵，$D_{ii}=d_i$，令$P=A{D}^{-1}$代表随机游走矩阵。$e_i$表示一个长度为n的向量，其中除了第i个元素其他都为0。e是全为1的向量。

导率：这个东西做社团发现的应该没有不知道的

矩阵的指数：这个东西在后面的推导中会用到：

3. FINDING SMALL CONDUCTANCE COMMUNITIES WITH DIFFUSIONS

图的扩散满足下面的性质：

这个式子表示从S出发进行随机游走，直到f在整个图上收敛。
其中${\sum }_k{\alpha }_k=1$，s是初始的种子向量。${\alpha }_k$递减保证扩散最终收敛。
值得注意的是，如果我们选取一个节点i作为种子，那么$s=e_i$，如果我们选取一团节点S作为种子节点，那么$s={\sum }_{i\in S}{e}_i/|S|$。

给定从一个初始向量s扩散而来的向量f，我们可以按照f中的值比上节点的度（惩罚那些度很大的节点）将节点排序（f中节点值越大说明越包含在由种子节点形成的社区中），我们最终可以使用这个序列找到能使得社区的导率最小的那些节点形成一个社区。

我们说明一下pagerank的扩散公式：

heat kernel其实只是将pagerank中的$\alpha$换成$t^k/k!$，那么公式就变成了

heat kernel与pagerank的比较：
heat kernel的系数$t^k/k!$收敛的比${\alpha }^k$更快，所以会将更大的权重放在距离种子节点更近的节点上。
下图可以看到heat kernel的收敛，确实比pagerank收敛地快很多。

4. algorithm

因为要准确地计算heat kernel的h向量需要将指数计算到无穷大，所以这是不可能的。我们希望近似计算heat kernel向量，不希望这个过程时间复杂度太高。

这里我们提出了一种名为hk-relax的近似算法：

回忆一下，我们寻找最优导率的社区的最后一步是将每个节点的h值除以他们的度，所以我们想要去计算一个近似向量$x\approx h$，满足下面的约束：

通过使用矩阵指数的标准性质，我们可以得到
所以我们的x向量满足
∣ ∣ D − 1 e − t e x p { t P } − D − 1 x ∣ ∣ < ϵ ||D^{-1}e^{-t}exp\{tP\}-D^{-1}x||<\epsilon ∣∣D−1e−texp{tP}−D−1x∣∣<ϵ
我们两边同时乘以$e^t$，得到我们希望找到一个向量y，对于所有的i，满足

4.1 Taylor Polynomial for exp{X}

我们首先使用泰勒多项式来近似估计矩阵的指数：
我们首先选择一个N使得

在这个前提下，如何我们能够找到y，满足：

那么根据三角形定理，我们有：

所以我们只需要找到近似估计$T_N(tP)s$的y就好了

4.3 Deriving a linear system

令$v_k$为$T_N(tP)s$的第k项，我们有：
$$T_N(tP)s=s+\frac{t}{1}Ps+\dots +\frac{t^N}{N!}{P}^{N}s$$

$$=v_0+v_1+\dots +v_N$$
我们注意到

令$v=[v_0;v_1;\dots v_N]$
一个近似的向量${\sum }_{k=0}^N{v}_k\approx T_n(tP)s$，也就是我们希望得到的$e^{t}h$的近似向量。

4.4 The hk-relax algorithm

有了上面的证明，我们可以这样计算：
我们保持一个最终的y向量，y向量开始是全0向量，我们再保持一个向量$v_k$，代表$T_N(tP)s$的第k项，$v_k$的初始值$v_0=s$，也就是是初始的种子向量，我们每次把$v_k$的值加入y中，再使用
来计算$v_{k+1}$，迭代这个过程直到y收敛。

我们基于这个思想对你这个过程又进行了一点改进，改进的原理是：
因为图中的节点数很多，$v_k$是一个稀疏向量，里面大多数值都是0，这些值对$v_{k+1}$的贡献值也是0，我们可以不用理会这些值。

所以我们可以这么做：
我们的数据结构如下：
y：我们要计算的diffusions向量
r:一个列表，列表的index是(u,k)，r[(v,k)]表示节点u在$v_k$中的值，也就是$v_k[u]$
q:一个队列，队列中保存的元素类型为(u,k)，也就是说保存那些$v_k[u]$不为0的值（我们对这里又做了一点改进，将所有$v_k[u]$小于某个阈值的元素舍去）。
最终的伪代码如下：

4.5 Choosing N

我们最后要解决的一个问题一就是如何选择N，在前面我们提到，我们要i选择一个N，使得

因为这个问题并不是这篇文章主要在注意的问题，我就直接把结论放上来：

直觉上可以看到N并不是很大，所以hk-relax的时间复杂度是可以保证的。

贴一段自己写的hk-relax

import networkx as nx
import numpy as np
import matplotlib.pyplot as plt


def hk_relax(g: nx.Graph, seeds, t, eps):
    N = get_N(t, eps)
    psis = get_psis(N, t)
    x = np.zeros(g.number_of_nodes())
    r = {}
    q = []
    for seed in seeds:
        r[(seed, 0)] = 1 / len(seeds)
        q.append((seed, 0))
    while len(q) > 0:
        (v, j) = q.pop(0)
        rvj = r[(v, j)]
        x[v] += rvj
        r[(v, j)] = 0
        mass = (t * rvj / (j + 1)) / g.degree[v]
        for u in g.neighbors(v):
            next = (u, j + 1)
            if j + 1 == N:
                x[u] += rvj / g.degree[v]
                continue
            if r.get(next, -1) == -1:
                r[next] = 0.
            thresh = np.exp(t) * eps * g.degree[u]
            thresh = thresh / (2 * N * psis[j + 1])
            if r[next] < thresh <= r[next] + mass:
                q.append(next)
            r[next] += mass
    return x


def get_N(t:float, eps):
    N = int(max(0,np.floor(t-2)+1))
    temp = t
    ans = eps
    while ans >= eps / 2:
        N += 1
        temp *= t / (N + 1)
        ans = temp * (N + 2) / (N + 2 - t)
    return N


def get_psis(N, t):
    temp = 1
    jiecheng = [temp]
    for i in range(1, N + 1):
        temp *= i
        jiecheng.append(temp)
    psis = []
    for j in range(N + 1):
        ans = 0.
        temp = 0.
        for m in range(N - j + 1):
            if m == 0:
                temp = 1
            else:
                temp *= t
            ans += (jiecheng[j] / jiecheng[m + j]) * temp
        psis.append(ans)
    return psis

g=nx.karate_club_graph()
x=hk_relax(g,[32],3,0.0001)
nx.draw(g,with_labels=True)
plt.show()

得到x之后，每个点的x值除以度排个序，从前向后选择能使得当前的节点组成的社区的导率最小的那些节点，就是我们得到的社区。