社群發現之標籤傳播演算法（LPA）python實現

郝hai發表於2024-04-26
演算法Python

社群發現在圖領域中備受關注，其根源可以追溯到子圖分割問題。在真實的社交網路中，使用者之間的聯絡緊密度不盡相同，導致形成了不同的社群結構。社群發現問題主要分為兩類：非重疊和重疊社群。非重疊社群發現指的是每個節點僅屬於一個社群，社群之間沒有交集。在非重疊社群發現中，有多種解決方法。其中，基於模組度的演算法透過最大化模組度來劃分社群，以找到最優的社群結構。另一種常見的方法是基於標籤傳播的演算法，例如標籤傳播演算法（LPA）。這種演算法透過在節點之間傳播標籤資訊，直到達到收斂條件，從而識別出社群結構。標籤傳播演算法相對簡單而且高效，適用於大規模圖的社群發現。這裡將重點討論標籤傳播演算法，探討其原理、優缺點以及在實際應用中的效果和侷限性。

一、標籤傳播演算法概述

標籤傳播演算法（Label Propagation Algorithm，簡稱LPA）是一種用於圖領域中的社群發現方法。其基本原理是利用節點之間標籤的傳播來將具有相似特徵的節點劃分到同一社群。演算法開始時，每個節點被賦予一個唯一的標籤，代表其所屬的潛在社群。隨後，每個節點根據其鄰居節點的標籤資訊進行標籤更新，迭代更新直至達到收斂狀態。具體而言，標籤傳播演算法的迭代過程如下：

初始化： 為每個節點分配一個唯一的標籤。
節點標籤更新： 隨機選擇一個節點，統計其鄰居節點的標籤分佈情況，選擇出現頻率最高的標籤作為該節點的新標籤。若有多個標籤出現頻率相同，則隨機選擇一個作為新標籤。
迭代： 重複節點標籤更新步驟，直到標籤分佈不再發生變化或達到預設的迭代次數。
透過這個過程，標籤傳播演算法能夠發現網路中的社群結構，將具有相似特徵的節點歸為同一社群。

1.1 標籤傳播演算法的優缺點

優點：
邏輯簡單： 演算法的核心思想直觀易懂，易於理解和實現。
無需預設社群數量： 不需要事先確定社群的數量，能夠根據網路結構自動劃分社群。
執行效率高： 演算法在迭代過程中只關注區域性標籤更新，避免了全域性最佳化問題，因此具有較高的執行效率。
缺點：
隨機性強： 演算法結果的穩定性可能受到初始化條件和迭代順序的影響。
無法處理重疊社群： 每個節點只能屬於一個社群，無法處理節點屬於多個社群的情況。
對噪聲和異常值敏感： 容易受到噪聲和異常值的影響，可能導致社群劃分不準確。

1.2標籤傳播演算法在實際應用中的效果與侷限性

標籤傳播演算法在實際應用中取得了一定的成功，尤其在處理大規模網路資料時表現出較高的效率。它可以幫助人們理解和分析複雜網路中的社交關係、使用者行為等。然而，演算法也存在一些侷限性：
穩定性問題： 演算法結果可能受到初始化和迭代順序的影響，導致結果不穩定。
無法處理重疊社群： 演算法無法準確處理節點同時屬於多個社群的情況，限制了其適用範圍。
對噪聲和異常值敏感： 在存在噪聲和異常值的情況下，演算法可能無法準確識別社群結構。
為了克服這些侷限性，研究者們提出了一些改進演算法，如COPRA、SLPA等。這些改進演算法透過引入更復雜的標籤傳播策略或額外的資訊來提高演算法的穩定性和準確性。
綜上標籤傳播演算法作為一種社群發現方法，具有簡單易懂、無需預設社群數量等優點，但也存在隨機性強、無法處理重疊社群等缺點。在實際應用中，需要根據具體情況選擇合適的演算法或進行改進以提高效果。隨著圖領域研究的不斷深入，我們期待看到更多針對標籤傳播演算法的改進和最佳化，使其在社群發現及其他相關領域發揮更大的作用。

二、標籤傳播演算法

2.1 演算法原理

LPA是基於圖的半監督方法，所以首先就需要對原資料進行圖方法的展示，如下圖：

那藍色的球（即節點）表示沒有標籤的樣本，棕色和黃色則可以表示有標籤的樣本。那樣本與樣本之間的連線（即權重）則可以表示它們之間的相似度。那麼這裡相似度的計算方法為：

\[w_{ij}=\exp^{\frac{\| \mathbf{x_i-x_j} \|^2}{\alpha^2}} \]

這裡的\(\alpha\)只是一個超引數而已。

當然這裡也只是說了這麼一種構建圖的方法，除此以外，也可以用KNN的方式來構建節點圖。不過KNN的方法會存在一種震盪的情況，導致資料的不穩定（主要是因為採用了同步更新法，換成非同步更新就沒有這個問題了）。不過KNN方法下的理論性沒有這種強，所以以這個方法展開來敘述LPA演算法。
有了上面的樣本與樣本間的相似度計算方法之後，我們可以構建一個基於相似矩陣的機率轉移矩陣\(P\)，\(P\)是一個N×N的矩陣：

\[P_{ij}=\frac{w_{ij}}{\sum_{k=1}^{n}w_{ik}} \]

顯然這裡是在做一個歸一化的行為。\(P_{ij}\)表示的就是節點\(i\)到\(j\)轉移的可能性，透過所有節點間轉移可能性的彙總，就可以得到我們需要的機率轉移矩陣P。
接著需要構建一個初始的狀態矩陣：

\[F = \begin{bmatrix} Y_l \\ Y_u \end{bmatrix} \]

上面\(Y_l\)是有標籤的樣本，下面的\(Y_u\)是沒有標籤的樣本。我們想要實現的就是把\(Y_l\)這些樣本也打上應有的標籤。既然有了機率轉移矩陣，也構建了初始的狀態矩陣，下面就是LPA演算法的實現過程了：

  • \(F=PF\)，透過這個計算，就可以把\(F\)的\(Y_u\)中部分無標籤資料打上了標籤。
  • 更新\(F\)中的\(Y_u\)樣本情況
  • 重複迭代上面過程這個過程，直到\(Y_u\)都有標籤為止。

思考一下：

既然我們關注點是\(Y_u\)部分，但是實際上我們每次在計算\(F=PF\)的過程中，\(Y_l\)都參與計算了。我們有沒有什麼更加簡便的計算方式，讓我們在計算過程中忽略\(Y_l\)呢？這是可以做到進一步最佳化：我們在構建機率轉移矩陣的時候，考慮把有標籤的樣本放在矩陣的前方，把無標籤的樣本放在矩陣的後方即可。那麼機率轉移矩陣就可以拆分為這個樣子了：

\[\begin{aligned} & P=\left[\begin{array}{ll} P_{L L} & P_{L U} \\ P_{U L} & P_{U U} \end{array}\right] \end{aligned} \]

這裡：
PLL表示已有標籤資料到其對應標籤的轉移機率，那肯定是1啊，所以PLL就是一個單位矩陣；
PLU表示已有標籤資料到其對不應標籤的轉移機率，那自然就是0了，所以PLU就是一個0矩陣；
PUL和PUU兩個矩陣則是未知的情況，因為沒有標籤的樣本究竟是屬於前面已經出現過的標籤，還是未知的標籤，目前還是未知狀態。
所以上面的P矩陣可以寫為：

\[\begin{aligned} P=\left[\begin{array}{cc} I & 0 \\ P_{U L} & P_{U U} \end{array}\right] \end{aligned} \]

繼續，我們採用這個\(P\)矩陣和\(F\)矩陣來實現上面\(LP\)演算法的迭代過程：

\[\begin{aligned} & \lim _{t \rightarrow \infty} P^t=\left[\begin{array}{cc} I & 0 \\ P_{U L} & P_{U U} \end{array}\right] \times\left[\begin{array}{cc} I & 0 \\ P_{U L} & P_{U U} \end{array}\right] \times \ldots \times\left[\begin{array}{cc} I & 0 \\ P_{U L} & P_{U U} \end{array}\right] \\ & =\left[\begin{array}{cc} I & 0 \\ \left(\sum_{t=0}^{\infty} P_{U U}^t\right) \cdot P_{U L} & P_{U U}^{\infty} \end{array}\right] \\ & =\left[\begin{array}{cc} I & 0 \\ \left(I-P_{U U}\right)^{-1} \cdot P_{U L} & 0 \end{array}\right] \\ & \end{aligned} \]

有兩點說明：

  • 正常情況下，未知的樣本不管屬於哪一個標籤，它的可能性是絕不會大於已知樣本屬於已知標籤的可能性（即1），所以：

\[P_{uu}<0\Rightarrow \lim_{t \rightarrow \infty}P_{\infty}^t=0 \]

  • 級數求和公式：

\[\sum_{k=0}^n x^k=\frac{1}{1-x} \quad|x|<1 \]

所以:

\[\left(\sum_{t=0}^{\infty} P_{U U}^t\right) \cdot P_{U L}=\left(I-P_{U U}\right)^{-1} \quad 0 \leq P_{U U}<1 \]

接著實現 \(F=P F\)，因為:

\[F=\left[\begin{array}{l} Y_L \\ Y_U \end{array}\right] \]

所以 \(F=P F\) 可以寫為:

\[\begin{bmatrix} \hat{Y}_L \\ \hat{Y}_U \end{bmatrix} = \lim_{t\rightarrow \infty}P^t =\begin{aligned} P=\left[\begin{array}{cc} I & 0 \\ (1-P_{UU})^{-1} \cdot P_{UL} & 0 \end{array}\right] \end{aligned} \begin{bmatrix} Y_L \\ Y_U \end{bmatrix} \]

最終得到:

\[\hat{Y}_U = (I - P_{UU})^{-1} \cdot P_{UL} \cdot Y_L \]

這樣，\(Y_L\)就可以完全排除我們的計算之外了。

所以對於LPA演算法的求解也有兩種辦法：

  • 正常求\(I-P_{UU}\)矩陣的逆，然後正常計算出結果即可。顯然這個在資料量較大時，比較耗費算力。
  • 迭代法，設定迭代的次數或者誤差的範圍，進行迭代即可。

2.2 演算法例項

再看個實際的例子吧：

接下來，按照上面的圖例來構建我們的機率轉移矩陣\(P\)：

\[P=\left[\begin{array}{ccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0.33 & 0 & 0.33 & 0 & 0.33 & 0 & 0 \\ 0 & 0 & 0 & 0.5 & 0 & 0.5 & 0 \\ 0 & 0.33 & 0 & 0 & 0.33 & 0 & 0.33 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{array}\right] \]

列表示的節點為\([1,7,2,3,4,5,6]\)，行表示為\([1,7,2,3,4,5,6]\)，如\(P\)的第一行：[1,0,0,0,0,0,0]表示：1號球屬於標籤1的可能性為1，屬於其餘標籤的可能性為0；P的第四行：[0.33,0,0.33,0,0.33,0,0]表示：3號球連線了1、2和4號三個球，故共有三種屬於的可能性，分別為0.33，所以在1、2以及4標籤處的連線可能性為0.33，其餘可能性為0，其他行的資料都是這麼計算而來。

接著，我們採用迭代法，計算\(P\)的無窮大次方：

import numpy as np
# 定義矩陣P
P = np.array([
    [1, 0, 0, 0, 0, 0, 0],
    [0, 1, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0],
    [0.33, 0, 0.33, 0, 0.33, 0, 0],
    [0, 0, 0, 0.5, 0, 0.5, 0],
    [0, 0.33, 0, 0, 0.33, 0, 0.33],
    [0, 0, 0, 0, 0, 1, 0]
], dtype=float)
# 計算矩陣P的n次冪
n=100
Pn = np.linalg.matrix_power(P, 100)
# 輸出結果
print(Pn)

計算的結果是：

\[P=\left[\begin{array}{ccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0.7316 & 0.239 & 0 & 0 & 0 & 0 & 0 \\ 0.7316 & 0.239 & 0 & 0 & 0 & 0 & 0 \\ 0.4853 & 0.4853 & 0 & 0 & 0 & 0 & 0 \\ 0.239 & 0.7316 & 0 & 0 & 0 & 0 & 0 \\ 0.239 & 0.7316 & 0 & 0 & 0 & 0 & 0 \end{array}\right] \]

再反過來看當時的\(P\)的求極限過程, 這個結果顯然驗證了當初的假設都是正確的。所以:

\[\begin{aligned} & \hat{Y}_U=\left(I-P_{U U}\right)^{-1} \cdot P_{U L} \times Y_L \\ & =\left[\begin{array}{cc} 0.7316 & 0.239 \\ 0.7316 & 0.239 \\ 0.4853 & 0.4853 \\ 0.239 & 0.7316 \\ 0.239 & 0.7316 \end{array}\right] \times\left[\begin{array}{ll} 1 & 0 \\ 0 & 1 \end{array}\right]=\left[\begin{array}{cc} 0.7316 & 0.239 \\ 0.7316 & 0.239 \\ 0.4853 & 0.4853 \\ 0.239 & 0.7316 \\ 0.239 & 0.7316 \end{array}\right] \end{aligned} \]

根據計算的結果:對於樣本1的橙色標籤來說，樣本2、3、4、5和6屬於它的可能性分別為：[0.7316,0.7316,0.4853,0.239,0.239]；對於樣本7的藍色標籤來說，樣本2、3、4、5和6屬於它的可能性分別為：[0.239,0.239,0.4853,0.7316,0.7316]，綜上：2和3是橙色；5和6是藍色；對於4來說，還是無法確定的。如圖：

除此以外，還有對一個對label propagation改進的一個演算法，即label spreading演算法。這個演算法本質上和label propagation是一樣的，不過它在迭代過程中，對label propagation的迭代方程進行了一個簡單的改進（就是加上了一個正則化），如公式：

\[F=P F \xrightarrow{\text { 正則化 }} F=\alpha P F+(1-\alpha) Y \]

這裡前後的\(P\)是有點區別的，但是基本一樣的方式計算得到的。不過重點不是這個機率轉移矩陣，而是折中引數α，它的取值在[0,1]之間，一般都是取0.2。

2.3 電商評價同步行為

  • 電商評價同步行為示例資料group

如第一行 x∩y 代表使用者 BUY_01 和 BUY_02 在同一個小時內對同一商戶進行評價的情況有 2 次，即他們出現同步行為的次數為 2 。

這裡同步行為次數是基於商戶（約束條件）去重作統計的，即A和B對同一個商家Y在不同時間段有多次匹配行為時，只算作一次。我們也可以考慮不對商家作去重處理，不過需要注意的是在傑斯卡相似度計算過程用到的\(x_cnt\)等指標邏輯也要保持一致，要麼都去重，要麼都不去重。如果採用不去重法得到的應是：

  • 將同步行為資料轉化成圖

由於 LPA 的演算法輸入為圖，需要將其轉換成圖形式

import matplotlib.pyplot as plt
import networkx as nx
tuples = [tuple(x) for x in group[['Buy_x','Buy_y','x∩y']].values] ##x∩y為weight
G = nx.Graph() ##本次採用無向圖，對應有DiGraph, MultiGraph, MultiDiGraph等畫圖選項
G.add_weighted_edges_from(tuples) ##如果是二元組，使用add_edges_from；weight另外定義
pos=nx.spring_layout(G)
nx.draw_networkx(G,pos)

示例樣本量太少，不易作社群劃分，手工新增多個社群節點

G.add_edge('BUY_04','BUY_05',weight=2)
G.add_edge('BUY_05','BUY_06',weight=1)
G.add_edge('BUY_05','BUY_07',weight=2)
G.add_edge('BUY_06','BUY_08',weight=1)
G.add_edge('BUY_07','BUY_08',weight=3)
G.add_edge('BUY_08','BUY_09',weight=1)
G.add_edge('BUY_08','BUY_10',weight=4)
G.add_edge('BUY_09','BUY_11',weight=2)
pos=nx.spring_layout(G)
labels = nx.get_edge_attributes(G,'weight')
nx.draw_networkx(G,pos)
nx.draw_networkx_edge_labels(G,pos,edge_labels=labels)

  • 使用 LPA 演算法劃分社群

使用同步行為出現的次數即 x∩y 值作為輸入權重

from networkx.algorithms.community import asyn_lpa_communities as lpa
# LPA本身不穩定，會存在波動
com = list(lpa(G,'weight')) #seed
print('社群數量',len(com))


  • 對劃分的社群進行業務打標

演算法將使用者分成了三個社群，可以將每個社群中的使用者的明細資料拉取出來，觀察他們的行為特徵。比如對於同一社群的使用者，觀察他們的評級時間密集程度如何，評價商戶有沒有集中性，評論相似性，所用裝置環境是否正常或者有共性等，基於業務考慮對社群進行打標，如人肉刷評群體。

三、標籤傳播演算法的Python算例



import numpy as np
import matplotlib.pyplot as plt
from sklearn.semi_supervised import LabelPropagation
import math

def show(Mat_Label, labels, Mat_Unlabel, unlabel_data_labels):
    for i in range(Mat_Label.shape[0]):
        if int(labels[i]) == 0:
            plt.plot(Mat_Label[i, 0], Mat_Label[i, 1], 'Dr')
        elif int(labels[i]) == 1:
            plt.plot(Mat_Label[i, 0], Mat_Label[i, 1], 'Db')
        else:
            plt.plot(Mat_Label[i, 0], Mat_Label[i, 1], 'Dy')

    for i in range(Mat_Unlabel.shape[0]):
        if int(unlabel_data_labels[i]) == 0:
            plt.plot(Mat_Unlabel[i, 0], Mat_Unlabel[i, 1], 'or')
        elif int(unlabel_data_labels[i]) == 1:
            plt.plot(Mat_Unlabel[i, 0], Mat_Unlabel[i, 1], 'ob')
        else:
            plt.plot(Mat_Unlabel[i, 0], Mat_Unlabel[i, 1], 'oy')

    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.xlim(0.0, 12.)
    plt.ylim(0.0, 12.)
    plt.show()

def loadCircleData(num_data):
    center = np.array([5.0 ,5.0])
    radiu_inner = 2
    radiu_outer = 4
    num_inner = int(num_data / 3)
    num_outer = num_data - num_inner
    data = []
    theta = 0.0
    for i in range(num_inner):
        pho = (theta % 360) * math.pi / 180
        tmp = np.zeros(2 ,np.float32)
        tmp[0] = radiu_inner * math.cos(pho) + np.random.rand(1) + center[0]
        tmp[1] = radiu_inner * math.sin(pho) + np.random.rand(1) + center[1]
        data.append(tmp)
        theta += 2
    theta = 0.0
    for i in range(num_outer):
        pho = (theta % 360) * math.pi / 180
        tmp = np.zeros(2, np.float32)
        tmp[0] = radiu_outer * math.cos(pho) + np.random.rand(1) + center[0]
        tmp[1] = radiu_outer * math.sin(pho) + np.random.rand(1) + center[1]
        data.append(tmp)
        theta += 1
    Mat_Label = np.zeros((2 ,2) ,np.float32)
    Mat_Label[0] = center + np.array([-radiu_inner + 0.5 ,0])
    Mat_Label[1] = center + np.array([-radiu_outer + 0.5 ,0])
    labels = [0 ,1]
    Mat_Unlabel = np.vstack(data)
    return Mat_Label ,labels ,Mat_Unlabel

def loadBandData(num_unlabel_samples):
    Mat_label = np.array([[5.0 ,2.] ,[5.0 ,8.0]])
    labels = [0 ,1]
    num_dim = Mat_label.shape[1]
    Mat_Unlabel = np.zeros((num_unlabel_samples ,num_dim) ,np.float32)
    Mat_Unlabel[:num_unlabel_samples // 2 ,:] = (np.random.rand(num_unlabel_samples // 2 ,num_dim) - 0.5) * np.array \
        ([3 ,1]) + Mat_label[0]
    Mat_Unlabel[num_unlabel_samples // 2:num_unlabel_samples ,:] = (np.random.rand(num_unlabel_samples // 2
                                                                                  ,num_dim) - 0.5) * np.array([3 ,1]) + Mat_label[1]
    return Mat_label ,labels ,Mat_Unlabel

if __name__ == "__main__":
    num_unlabel_samples = 2000  # 樣本數量增加以滿足最近鄰居的要求
    Mat_label, labels, Mat_Unlabel = loadCircleData(num_unlabel_samples)

    # 使用LabelPropagation演算法進行標籤傳播
    label_prop_model = LabelPropagation(kernel='knn', n_neighbors=min(10, num_unlabel_samples -1), max_iter=300)
    label_prop_model.fit(np.vstack((Mat_label, Mat_Unlabel)), labels + [-1] * num_unlabel_samples)
    unlabel_data_labels = label_prop_model.transduction_

    # 顯示結果
    show(Mat_label, labels, Mat_Unlabel, unlabel_data_labels)

總結

標籤傳播演算法（Label Propagation Algorithm，LPA）是一種用於社群發現的簡單而有效的半監督方法，其核心思想是透過節點之間標籤的傳播來劃分社群，將具有相似特徵的節點歸為同一社群。演算法適用於各種型別的圖，包括社交網路、通訊網路等。在實際應用中，標籤傳播演算法被廣泛應用於社交網路分析、推薦系統、生物資訊學等領域。例如，在社交網路中，它可以幫助識別使用者群體、發現潛在的社交圈子；在推薦系統中，可以利用使用者行為資料進行社群劃分，從而提高推薦的準確性和個性化程度；在生物資訊學中，可以用於基因調控網路的模組化分析，發現基因間的關聯性。此外，標籤傳播演算法還可用於影像分割、文字聚類等領域。例如，在影像分割中，可以將畫素點作為圖的節點，利用它們的相似性進行社群劃分，從而實現影像的自動分割；在文字聚類中，可以將文件表示為圖的節點，透過詞語之間的關聯性進行社群劃分，實現文字的自動分類和聚類。總之標籤傳播演算法是一種簡單而有效的社群發現方法，在各種領域都有著廣泛的應用前景。透過利用其優勢並克服其侷限性，可以為社交網路分析、推薦系統、生物資訊學等領域的問題提供有力的解決方案。

參考文獻

  1. 兩種 Python 方法實現社群發現之標籤傳播演算法（LPA）
  2. LPA社群發現演算法介紹和程式碼示例
  3. 27.標籤傳播演算法

相關文章