前回に引き続き、Louvain法の話です。前回はPythonで使う方法を紹介しましたが、今回はこの手法が利用しているモジュラリティ(modularity)という指標について検証していきます。そもそもなぜこの式でグラフのノードのクラスタを評価できるのかってのを確認しましょう。
(ちなみに、モジュラリティを最大化するアルゴリズムについては今回の記事では取り上げません。あくまでもモジュラリティの定義について見ていきます。)
前回の記事でも書いてるのですが今回の記事では主役なので定義式を再掲します。
$$Q=\frac{1}{2m}\sum_{ij}\left[A_{ij}\ -\ \frac{k_ik_j}{2m} \right]\delta(c_i, c_j).$$
$A_{ij}$: ノード$i$, ノード$j$間のエッジの重み。
$k_i, k_j$: ノード$i$, ノード$j$それぞれに接続されたエッジの重みの合計。
$m$: グラフの全てのエッジの重みの合計。
$c_i, c_j$: ノード$i$, ノード$j$それぞれが所属しているコミュニティ。
$\delta$: クロネッカーのデルタ。二つの引数が等しければ1でそれ以外は0を返す関数。
一見しただけだと、これでクラスタを評価できるってわかりにくいですよね。
順番に見ていきましょう。
まず、$\delta(c_i, c_j)$の部分です。これ、$c_i$と$c_j$が別のクラスタだったら値が$0$なので、$\sum$のそれらの項は消滅し、さらに同じクラスタだったら値が$1$になるので、それらの項はそのままの値で残ります。
なのでこの式はこうなります。
$$Q=\frac{1}{2m}\sum_{c_iとc_jが同じクラスタとなるi,j}\left[A_{ij}\ -\ \frac{k_ik_j}{2m} \right].$$
そして、$A_{ij}$はエッジの重みなので、和の中身の前半部分だけに着目すると、次の様に言えます。
$$\frac{1}{2m}\sum_{c_iとc_jが同じクラスタとなるi,j}A_{ij} = \frac{1}{2m}\{クラスタ内部のエッジの重みの総和\}$$
$m$はクラスタに関係なく、グラフによって定まってる定数なので無視して良いですから、この部分は確かに同じクラスタの中にたくさんエッジがあると値が大きくなり評価指標としての意味を持っている様に見えます。
ただし、この項しかないと、全部のノードを1個のクラスタとした場合に$Q$が最大となっしまうので、それに対応する必要があります。そこで登場するのが残りの項です。
もし、ノード$i$とノード$j$の間にエッジがなかったりあってもウェイトが非常に小さかったりすると、負の数を足すことによってモジュラリティが下がる様になってるのです。たとえば、エッジがない場合、$A_{ij}$が0ですから、$A_{ij}\ -\ \frac{k_ik_j}{2m} = -\frac{k_ik_j}{2m} $です。
そして、$k_i$、$k_j$はそれぞれのノードに接続されたエッジの重みの和ですから、重みの大きなエッジをたくさん持っているノード間についてはこの値の絶対値は大きくなります。
つまり、エッジがたくさんあるノード同士にも関わらずそれらのノード間にエッジがなかったら強くペナルティを課すよ、ってのがこの項の意味なのです。
分母の$2m$とかにもちゃんと意味があるのですが細かい話になるので割愛します。これらの細かい工夫により、-0.5から1の値を取る指標となっています。計算量とのバランスもとりながら非常によく考えられた指標だと思いました。
最後に、前回の記事でも少しだけ触れているのですが、この$\sum_{ij}$が取る和の$ij$部分について確認したのでその話書いておきます。
これは$i$と$j$がそれぞれ独立に全てのパターンを取ります。$(i, j)$と$(j, i)$は個別に足されるし、$(i, i)$も考慮されるってことですね。ライブラリの結果と、自分で実装した結果を並べてみて確認しました。とりあえず復習兼ねて前回のサンプルのグラフでやります。まずライブラリ利用。
import numpy as np
import networkx as nx
import community
# 前回の記事で生成したのと同じグラフを再現する。
# 30個のノード
node_list = list(range(30))
# 同じクラスタ内は0.5, 別クラスタは0.02の確率でエッジを生成する
edge_list = []
np.random.seed(1) # シード固定
for i in range(30):
for j in range(i+1, 30):
if i // 10 == j // 10:
if np.random.rand() < 0.5:
edge_list.append((i, j))
else:
if np.random.rand() < 0.02:
edge_list.append((i, j))
# グラフ生成
G = nx.Graph()
G.add_nodes_from(node_list)
G.add_edges_from(edge_list)
# コミュニティーの検出
partition = community.best_partition(G)
# ライブラリによるモジュラリティの計算
print(community.modularity(partition, G))
# 0.5316751700680272同じ値を自分で計算して出してみましょう。エッジのウェイトが全部1なので、$m$はただのエッジの本数になるし、$k_i$, $k_j$はそれぞれのノードの次数に等しいことに注意してください。
Q = 0
m = G.size()
# iと jは全部の組み合わせをとる。
for i in range(len(G)):
for j in range(len(G)):
if partition[i] == partition[j]:
# ノードi と ノードjの間にエッジがあったら1(ウェイト)を足す
if (i, j) in G.edges:
Q += 1
# Σの後半部分
Q -= G.degree(i) * G.degree(j) / (2*m)
# 全体を2mで割る
Q/=2*m
print(Q)
# 0.5316751700680273一致してますね。
ここから余談ですが、Qは -0.5から 1の値を取りますが、どんなグラフについてもQが-0.5を取ったり1を取ったりする分割が存在するわけではなさそうです。上記のサンプルでも0.53程度に留まっていますしね。
-0.5 の方は、完全２部グラフを構築して、さらにその2部をそれぞれコミュニティとすることで、同じクラスタ内部には1本もエッジがないのに別クラスタのノード間には必ずエッジが存在するという状態にすると発生します。
また、逆にQが大きくなるのは非連結な多数の完全グラフからなるグラフをそれぞれの完全グラフを一つのクラスタとすると、非連結なグラフが増えるにつれて1に近づいていくのを確認できています。ただ、近づくだけで1になることはないんじゃないでしょうか。(証明はできていませんが。)
便利なのでなんとなく使っていた手法でしたが、改めて指標を検証し複数のグラフや分割について自分で計算してみたことで理解を深めることができました。
モジュラリティは定義中に記号多いですしぱっと見難しそうに見えるのですが、落ち着いて見れば四則演算だけで構成されている非常にシンプルな指標なので、興味があれば皆さんもいろいろ試して遊んでみてください。