#概要
コミュニティのクラスタリングがいずれ仕事にも使えそうなので、お試ししてみた備忘録を残す。
いろいろ調べててGraphRicciCurvatureを見つけた。ユークリッド空間ではなく位相空間で引っ張ったり伸ばしたりするやつらしい。なんか面白そう。
でもGraphRicciCurvatureで検索してもヒットしなかったので、Qitta初か？
それほど新しいアルゴリズムではないようだが、日本語情報はいつにも増して少ない。
- 実施期間: 2021年8月
- 環境：Ubuntu20.04 LTS
- パケージ：GraphRicciCurvature, networkxなど
##1. 準備
パケージをインストールする。前回のpostが動作するパケージはインストール済みとする。
オフィシャルの通りにインストールするだけだが、依存パケージのコンパイルの必要があるので諸々必要。
# 入っていなければ。
sudo apt install g++
sudo apt install cmake
sudo apt install clang
# 各パケージのオフィシャルの通り
pip3 install networkit
pip install POT
pip3 install GraphRicciCurvature
今回使用するパケージをすべてimportする。
今回もお試しのデータセットはkarateを使用する。
基本的にはオフィシャルチュートリアルに沿ってかいつまんで記載する。
チュートリアルはColabでも動作するように書かれているが、何故かcommunityが見つからないとエラーが出る。Colab上でGithubのソースからpython-louvainをインストールすれば良いような気もするがGPUは必要ないので試していない。
import networkx as nx
import numpy as np
import math
# matplotlib setting
%matplotlib inline
import matplotlib.pyplot as plt
# to print logs in jupyter notebook
import logging
logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.ERROR)
# load GraphRicciCuravture package
from GraphRicciCurvature.OllivierRicci import OllivierRicci
from GraphRicciCurvature.FormanRicci import FormanRicci
# load python-louvain for modularity computation
import community as community_louvain
# for ARI computation
from sklearn import preprocessing, metrics
import GraphRicciCurvature
print(GraphRicciCurvature.__version__)
0.5.3
helper関数は長いので最後に記載する。
##2. GraphRicciCurvatureでクラスタリング
トポロジーの知識が必要なようだが、そんなものは専攻したことがないので理屈は理解できない。
パケージの解説によると、
「edgeのRicci曲率がグラフの構造に重要な役割を果たしていることが観察されています。正の曲率を持つedgeはクラスタ内のedgeを表し、負の曲率を持つedgeはクラスタ内の架け橋となります。また、負の曲率を持つedgeはグラフの連結性に大きく関係しており、連結されたグラフから負の曲率を持つedge取り除くと、グラフはすぐに切断されてしまいます。
Ricci flowは、グラフのエッジのRicci曲率を均一化するプロセスである。Ricci flowは、与えられたグラフに対して、各edgeに"Ricci flow metric"をedge weightとして与え、このedge weightより小さいければグラフのRicci曲率がどこでもほぼ等しくなるようにするものです。この"Ricci flow metric"がコミュニティを検出できます。
Ricci曲率とRicci flow metricの両方とも、グラフ分類のためのグラフフィンガープリントとして機能します。グラフが異なれば、edgeのRicci曲率分布もRicci flow metricも異なります。」
あわせてココのレクチャによると、metric tensor: gは球体のどの部分を伸長・収縮させるかを示す。その曲率をRicci curveture: Rと呼び、一般式はココに定義あり。
従い、gがわかればRもわかる。
これをRicci flowに適用する。Ricci flowではgをiterativeに変化させ球体を球に近づける。
このとき球体の凹んだ部分(R<0)は膨らみ、凸った部分(R>0)はへこむ。つまり凹んだ部分では膨らますためgは増加し、凸った部分ではgは減少する。
R = -dg/dt
従いほっておくとR>0なら縮退し消滅する。
縮退の過程で球体の一部が限りなくゼロに近づく（ネック）が放置するとゼロに縮退するのでこのネック部分をマニュアルで削除し、その部分に球体をかぶせる(Surgery Theory)。この後の解説がどうしても納得できなかった。
まぁ、このパケージではこのRがiterationを通じて一定になるようにしているのでグラフは消滅しないし、更新されるedge weightを閾値としてクラスタ間を切断(surgery)することもできる。
アプローチが異なる2つの手法が用意されている。
Ollivier-Ricci curvature:
最適輸送理論に基づく曲率。各edgeの曲率をミクロに見ることができますが、計算量が多くなります。
Forman-Ricci curvature:
CWコンプレックスに基づく曲率です。各エッジの曲率をマクロ的に見ることができ、計算速度も非常に軽い。
使い方は同じなので前者のコードを残す。
全APIはココを参照
###2.1 Ollivier-Ricci curvature
G = nx.karate_club_graph()
draw_graph(G_rf)
orc = OllivierRicci(G, alpha=0.5, verbose="INFO")
orc.compute_ricci_curvature()
G_orc = orc.G.copy() # save an intermediate result
show_results(G_orc)
Karate Club Graph, first 5 edges:
Ricci curvature of edge (0,1) is 0.111111
Ricci curvature of edge (0,2) is -0.143750
Ricci curvature of edge (0,3) is 0.041667
Ricci curvature of edge (0,4) is -0.114583
Ricci curvature of edge (0,5) is -0.281250
通常、Ricci curvature(以降R)は[-1, 1]の間に入る。
weightは重要なパラメータだが、元のに設定されていなければdefaultで1となる。
###2.2 Ricci flow
Rを一定に、かつ曲率をなめらかにしながらedge weightを更新する。Rが大きくマイナス/プラスになるとそのedgeを伸ばす/縮めることとなる。
経験的には20~50回iterationすると、Rはなめらかになるらしい。
とりあえず一度iterationを回すとどうなるか見てみる。
# Start a Ricci flow with Lin-Yau's probability distribution setting with 4 process.
orf = OllivierRicci(G, alpha=0.5, base=1, exp_power=0, proc=4, verbose="INFO")
# Do Ricci flow for 2 iterations
orf.compute_ricci_flow(iterations=2)
orf.set_verbose("ERROR") # mute logs
orf.compute_ricci_flow(iterations=50) # 上の2回(iterations=2)に続き、行なうことになる。
G_rf = orf.G.copy()
show_results(G_rf)
Karate Club Graph, first 5 edges:
Ricci curvature of edge (0,1) is -0.002740
Ricci curvature of edge (0,2) is -0.002740
Ricci curvature of edge (0,3) is -0.002740
Ricci curvature of edge (0,4) is -0.002740
Ricci curvature of edge (0,5) is -0.002740
Rが一定値に収束し、全edge weightが調整されていることがわかる。
draw_graph(G_rf)
surgery実行していないので全edgeはそのまま。色は元の"club"属性を示す（以降同じ）。
###2.3 Ricci flow for Community detection
クラスタリングにおいては、weight閾値を見つけそれよりも大きなweightを持つedgeを取り除く。
例えば,weight=1.5と1.0でカットした場合を見てみる。
draw_graph(my_surgery(G_rf, cut=1.5))
左上のnode=[4,5,6,10,15]は前回のLouvainのデンドロが最初に見つけた小グループと同じである。全く違うアルゴリズムだが見つけるものは見つけてくれる。
draw_graph(my_surgery(G_rf, cut=1.0))
cut=1.5や1.0はaccuracyを見ながら最適値に調整する。次項で説明する。
###2.4 Optimal transportation and Communities
別のインスタンス(orf2)を作成して、それでやってみる。
draw_graph(my_surgery(G_rf, cut=1.0))
orf2 = OllivierRicci(G, alpha=0.5, base=math.e, exp_power=1, verbose="ERROR")
orf2.compute_ricci_flow(iterations=50)
G_rf2 = orf2.G.copy()
show_results(G_rf2)
Karate Club Graph, first 5 edges:
Ricci curvature of edge (0,1) is -0.019349
Ricci curvature of edge (0,2) is -0.001516
Ricci curvature of edge (0,3) is -0.010322
Ricci curvature of edge (0,4) is -0.001809
Ricci curvature of edge (0,5) is -0.002310
ヒストグラムが前項と異なるのは引数がbase=math.e, exp_power=1に変わっているから。これとalphaを調整する必要がありそう。チュートリアルによると、もしくはbase=math.e, exp_power=2が理論的に最適（熱拡散に近い）になるとのこと。なお、exp_power=0にすると元論文となる。
alphaは、確率分布の元nodeと隣接nodeの分割割合になる。[0, 1.0]
このときのaccuracyを確認する。
check_accuracy(G_rf2, clustering_label="club",plot_cut=True)
チャート中のGood cutoffは多分weight=0から増加させながらaccが最初に落ちたときを拾っているぽい。
全体を見るとaccが最高なのは3.0 ~ 3.88の間なので3.88でcutしてみる。
draw_graph(my_surgery(G_rf2, cut=3.88))
うまく分離できている。node=9とnode=8を取り違えているのは前回のLouvainと同じだが、大きく異なるのはRによりedgeを任意に取り除くことができる点で、こっちのほうがよい。
ただ、3.0~3.88を求めるロジックやalphaの調整ロジックは考えないといけない。
一意に決まらないようであれば使うのは難しいか？
ちなみにcutoffを2.0にすると例の小クラスタが分離される。
draw_graph(my_surgery(G_rf2, cut=2.0))
答えはないのでクラスタリングを何に使うかによって調整しなければならない。
2.5 Helper関数
helper関数をオフィシャルから転記しておく。
def check_accuracy(G_origin, weight="weight", clustering_label="value", plot_cut=True):
"""To check the clustering quality while cut the edges with weight using different threshold
Parameters
----------
G_origin : NetworkX graph
A graph with ``weight`` as Ricci flow metric to cut.
weight: float
The edge weight used as Ricci flow metric. (Default value = "weight")
clustering_label : str
Node attribute name for ground truth.
plot_cut: bool
To plot the good guessed cut or not.
"""
G = G_origin.copy()
modularity, ari = [], []
maxw = max(nx.get_edge_attributes(G, weight).values())
cutoff_range = np.arange(maxw, 1, -0.025)
print(maxw)
for cutoff in cutoff_range:
edge_trim_list = []
for n1, n2 in G.edges():
if G[n1][n2][weight] > cutoff:
edge_trim_list.append((n1, n2))
G.remove_edges_from(edge_trim_list)
# Get connected component after cut as clustering
clustering = {c: idx for idx, comp in enumerate(nx.connected_components(G)) for c in comp}
# Compute modularity and ari
modularity.append(community_louvain.modularity(clustering, G, weight))
ari.append(ARI(G, clustering, clustering_label=clustering_label))
plt.xlim(maxw, 0)
plt.xlabel("Edge weight cutoff")
plt.plot(cutoff_range, modularity, alpha=0.8)
plt.plot(cutoff_range, ari, alpha=0.8)
if plot_cut:
good_cut = -1
mod_last = modularity[-1]
drop_threshold = 0.01 # at least drop this much to considered as a drop for good_cut
# check drop from 1 -> maxw
for i in range(len(modularity) - 1, 0, -1):
mod_now = modularity[i]
if mod_last > mod_now > 1e-4 and abs(mod_last - mod_now) / mod_last > drop_threshold:
if good_cut != -1:
print("Other cut:%f, diff:%f, mod_now:%f, mod_last:%f, ari:%f" % (
cutoff_range[i + 1], mod_last - mod_now, mod_now, mod_last, ari[i + 1]))
else:
good_cut = cutoff_range[i + 1]
print("*Good Cut:%f, diff:%f, mod_now:%f, mod_last:%f, ari:%f" % (
good_cut, mod_last - mod_now, mod_now, mod_last, ari[i + 1]))
mod_last = mod_now
plt.axvline(x=good_cut, color="red")
plt.legend(['Modularity', 'Adjust Rand Index', 'Good cut'])
else:
plt.legend(['Modularity', 'Adjust Rand Index'])
def show_results(G, curvature="ricciCurvature"):
# Print the first five results
print("Karate Club Graph, first 5 edges: ")
for n1,n2 in list(G.edges())[:5]:
print("Ricci curvature of edge (%s,%s) is %f" % (n1 ,n2, G[n1][n2][curvature]))
# Plot the histogram of Ricci curvatures
plt.subplot(2, 1, 1)
ricci_curvtures = nx.get_edge_attributes(G, curvature).values()
plt.hist(ricci_curvtures,bins=20)
plt.xlabel('Ricci curvature')
plt.title("Histogram of Ricci Curvatures (Karate Club)")
# Plot the histogram of edge weights
plt.subplot(2, 1, 2)
weights = nx.get_edge_attributes(G, "weight").values()
plt.hist(weights,bins=20)
plt.xlabel('Edge weight')
plt.title("Histogram of Edge weights (Karate Club)")
plt.tight_layout()
def ARI(G, clustering, clustering_label="club"):
"""
Computer the Adjust Rand Index (clustering accuracy) of "clustering" with "clustering_label" as ground truth.
Parameters
----------
G : NetworkX graph
A given NetworkX graph with node attribute "clustering_label" as ground truth.
clustering : dict or list or list of set
Predicted community clustering.
clustering_label : str
Node attribute name for ground truth.
Returns
-------
ari : float
Adjust Rand Index for predicted community.
"""
complex_list = nx.get_node_attributes(G, clustering_label)
le = preprocessing.LabelEncoder()
y_true = le.fit_transform(list(complex_list.values()))
if isinstance(clustering, dict):
# python-louvain partition format
y_pred = np.array([clustering[v] for v in complex_list.keys()])
elif isinstance(clustering[0], set):
# networkx partition format
predict_dict = {c: idx for idx, comp in enumerate(clustering) for c in comp}
y_pred = np.array([predict_dict[v] for v in complex_list.keys()])
elif isinstance(clustering, list):
# sklearn partition format
y_pred = clustering
else:
return -1
return metrics.adjusted_rand_score(y_true, y_pred)
def my_surgery(G_origin: nx.Graph(), weight="weight", cut=0):
"""A simple surgery function that remove the edges with weight above a threshold
Parameters
----------
G_origin : NetworkX graph
A graph with ``weight`` as Ricci flow metric to cut.
weight: str
The edge weight used as Ricci flow metric. (Default value = "weight")
cut: float
Manually assigned cutoff point.
Returns
-------
G : NetworkX graph
A graph after surgery.
"""
G = G_origin.copy()
w = nx.get_edge_attributes(G, weight)
assert cut >= 0, "Cut value should be greater than 0."
if not cut:
cut = (max(w.values()) - 1.0) * 0.6 + 1.0 # Guess a cut point as default
to_cut = []
for n1, n2 in G.edges():
if G[n1][n2][weight] > cut:
to_cut.append((n1, n2))
print("*************** Surgery time ****************")
print("* Cut %d edges." % len(to_cut))
G.remove_edges_from(to_cut)
print("* Number of nodes now: %d" % G.number_of_nodes())
print("* Number of edges now: %d" % G.number_of_edges())
cc = list(nx.connected_components(G))
print("* Modularity now: %f " % nx.algorithms.community.quality.modularity(G, cc))
print("* ARI now: %f " % ARI(G, cc))
print("*********************************************")
return G
def draw_graph(G, clustering_label="club"):
"""
A helper function to draw a nx graph with community.
"""
complex_list = nx.get_node_attributes(G, clustering_label)
le = preprocessing.LabelEncoder()
node_color = le.fit_transform(list(complex_list.values()))
nx.draw_spring(G, nodelist=G.nodes(),
node_color=node_color,
cmap=plt.cm.rainbow,
alpha=0.8, with_labels=True)
@article{ni2019community,
title={Community detection on networks with ricci flow},
author={Ni, Chien-Chun and Lin, Yu-Yao and Luo, Feng and Gao, Jie},
journal={Scientific reports},
volume={9},
number={1},
pages={1--12},
year={2019},
publisher={Nature Publishing Group}
}
以上