1.はじめに
この記事は， pythonによるネットワーク分析（コミュニティ抽出），ネットワーク分析に便利なライブラリである NetworkXの簡単な使い方のまとめです．
コミュニティ抽出のアルゴリズムとして， スペクトラルクラスタリングをnumpyで実装します．
個人的なモチベーションとしては，この記事（Probabilistic Matrix Factorization を導出して Edward で実装する）を読んで，この本（関係データ学習）面白そうだなと思って読み始めたのでそのまとめのようなものです．
2.環境構築
2.1.今回の実施環境
- Python 3.6.0 :: Anaconda 4.3.1 (x86_64)
- numpy 1.11.3
2.2.インストール
pip install networkx==1.9.1
networkxの1.9.1をインストールします．最新版は1.11ですが，それを使うとデータをロードするときにNetworkXError: cannot tokenize 'graph' at (2, 1)が出ました．調べてみると「俺もそれ出たけど1.9.1にダウングレードしたら使えるようになったぜ」という人がいたので，それに習います．
3.データについて
3.1.データのソース
ここにネットワーク分析関係のお試しデータがいくつも公開されているので，この中から，
Les Miserables: coappearance network of characters in the novel Les Miserables. Please cite D. E. Knuth, The Stanford GraphBase: A Platform for Combinatorial Computing, Addison-Wesley, Reading, MA (1993).
を使います．
3.2.データの内容
このデータは，レミゼラブルの小説の中で，同じシーンに一緒に現れた人のデータだそうです．一緒に現れた回数がエッジの値になり，ノードは登場人物になります．
面白いデータだなあ，こんなデータあったのかあと思っていたら，コードを書き終わった後に，twitter界隈で有名なTJO氏のこの記事を見つけました．
UCI機械学習リポジトリのデータ（など）で遊ぶ(2)：『レ・ミゼラブル』の人物相関図
同じコミュニティ抽出の課題を，別の手法で試しています．このレミゼデータがどんなデータで，コミュニティ抽出をするとどんなことがわかるのかというのはこの記事を読むと大体わかるのではと思います．
この記事を読んだ感じ，ネットワーク分析のライブラリはpythonよりもRの方が揃っているのではと思いました．もしくはcytoscapeなどを使うべきなのでしょうか．（もしくは僕が他のものを知らないだけ）
3.3.データの形式
.gmlという拡張子のデータが，上記リンクのzipの中に入っています．gmlとは，Graph Modelling Languageの略で，グラフ構造のデータの保存形式です．
以下のような構造で，グラフのノードとエッジの情報を格納しています．
$ cat lesmis.gml
Creator "Mark Newman on Fri Jul 21 12:44:53 2006"
graph
[
node
[
id 0
label "Myriel"
]
node
[
id 1
label "Napoleon"
]
...
edge
[
source 14
target 11
value 1
]
edge
[
source 15
target 11
value 1
3.4.データの読み込み
.gmlデータは，NetworkXのread_gml関数でNetworkXの定義するグラフオブジェクトに読み込めます．
import networkx as nx
data_path = "./data/lesmis/lesmis.gml"
G = nx.read_gml(data_path)
これをnumpyの配列にしたい時は，
X = np.array(nx.to_numpy_matrix(G))
というようにできます．to_numpy_matrixの返り値はnumpy matrixですが，個人的にndarrayの方が使いやすいので上のようにndarrayに直しています．
ちなみに，今回のスペクトラルクラスタリングでは，グラフの行列が対称行列であることを仮定しています．データをちょっと見ればわかることですが，NetworkXでは
nx.is_directed(G)
# False
という関数で，グラフが有向グラフ（対称行列でない）か無向グラフ（対称行列）かチェックできます．有向かどうかチェックしてFalseなので，無向グラフです．
4.コード
4.1.概要
ここでは，
- 非正規化グラフラプラシアンに基づくスペクトラルクラスタリング
- 酔歩正規化グラフラプラシアンに基づくスペクトラルクラスタリング
の二つを実装しています．
非正規化グラフラプラシアンに基づく方法では，クラスタリングを実行すると，ある一つのノードとそれ以外というようにクラス分けする傾向があり，クラスタリングとして微妙なので，その改良版が正規化グラフラプラシアンに基づくものです．
正規化グラフラプラシアンは正規化の仕方によって種類があり，その中の酔歩正規化グラフラプラシアンを実装しています．
4.2.クラスや関数の定義
実際にスペクトラルクラスタリングを行うクラスと，実行結果をプロットするのに必要な関数を定義します．
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from sklearn.cluster import KMeans
class graph_clustering():
def __init__(self, K, n_clusters, lap_kind):
"""
K : int
n_clusters : int
lap_kind : None or str
Ex. "rw"
"""
self.lap_kind = lap_kind
self.K = K
self.n_clusters = n_clusters
def _calc_g_laplacian(self, X):
D = np.diag(np.sum(X, axis=0))
L = D - X
self.L = L
self.D = D
return L
def _calc_rw_g_laplacian(self, X):
L = self._calc_g_laplacian(X)
Lrw = np.dot(np.linalg.inv(self.D), L)
self.L = Lrw
return Lrw
def _get_K_eigen_vec(self, Lam, V):
sorted_index = np.argsort(Lam.real)
Lam_K = Lam[sorted_index][0:self.K]
V_K = V[sorted_index][0:self.K]
self.Lam_K = Lam_K
self.V_K = V_K
return Lam_K, V_K
def _Kmeans_V_K(self, V_K):
kmeans = KMeans(n_clusters=self.n_clusters, random_state=0).fit(V_K.T.real)
clusters = kmeans.labels_
self.clusters = clusters
return clusters
def laplacian_clustering(self, X):
"""
X : ndarray
a matrix representation of undirected graph
"""
if self.lap_kind is None:
L = self._calc_g_laplacian(X)
elif self.lap_kind == "rw":
L = self._calc_rw_g_laplacian(X)
else:
raise NotImplementedError
Lam, V = np.linalg.eig(L)
Lam_K, V_K = self._get_K_eigen_vec(Lam, V)
clusters = self._Kmeans_V_K(V_K)
return clusters
# for plot
def get_labels_dic(G):
"""
G : networkx.classes.graph.Graph
"""
labels_dic = { key:G.node[key]['label'] for key in G.node.keys() }
return labels_dic
def get_labels_list(G):
labels_list = np.array([ G.node[key]['label'] for key in G.node.keys() ])
return labels_list
def split_labels(labels_list, clusters):
"""
labels_list : list
return value of get_labels_list function.
clusters : ndarray
return value of graph_clustering.laplacian_clustering
"""
class_uniq = np.unique(clusters)
node_index_split= []
labels_split = []
index = np.arange(clusters.shape[0])
for class_elem in class_uniq:
node_index_split.append(list(index[clusters == class_elem]))
labels_split.append(list(labels_list[clusters == class_elem]))
return node_index_split, labels_split
def plot_clusters(G, node_index_split):
"""
G : networkx.classes.graph.Graph
node_index_split : list (2-dim list)
return value of split_labels function.
"""
labels_dic = get_labels_dic(G)
pos = nx.spring_layout(G)
colors = [ cm.jet(x) for x in np.arange(0, 1, 1/len(node_index_split)) ]
for i, nodes in enumerate(node_index_split):
nx.draw_networkx(G, node_color=colors[i], labels=labels_dic, nodelist=nodes, pos=pos, font_color="r")
plt.show()
コードの詳細な説明は後回しにして，次のようにしてスペクトラルクラスタリングとその結果のプロットを行えます．
4.3.コードを動かす
data_path = "./data/lesmis/lesmis.gml"
G = nx.read_gml(data_path)
X = np.array(nx.to_numpy_matrix(G))
GClus = graph_clustering(K=15, n_clusters=5, lap_kind="rw")
clusters = GClus.laplacian_clustering(X)
labels_list = get_labels_list(G)
node_index_split, labels_split = split_labels(labels_list, clusters)
plot_clusters(G, node_index_split)
という感じで使います．
graph_clustering(K=15, n_clusters=5, lap_kind="rw")の引数は，Kが固有ベクトルを何個取るか，n_clustersは何個のクラスターに分割するか，lap_kindはスペクトラルクラスタリングのアルゴリズムの種類を与えます．
lap_kind=Noneで非正規化，lap_kind="rw"で酔歩正規化グラフラプラシアンのクラスタリングを実行します．
4.4.プロットするコードについて
少し工夫が必要であったところが，プロット部分です．plot_clusters()内のnx.draw_networkx関数によってプロットしていますが，識別されたクラスターごとに色を変えてプロットできるようなオプションがなかったため，このようにプロットしています．
nx.draw_networkx(G, node_color=colors[i], labels=labels_dic, nodelist=nodes, pos=pos, font_color="r")
このnodelistという引数でグラフ全体の中から部分的にプロットできたので，ノードをクラスターごとに分割し，node_color=colors[i]で色を変えつつプロットしています．
また，これはプロットをどんどん重ね書きするので，pos引数でポジションを指定してあげないと，毎回ポジションが変わってしまってぐちゃぐちゃになってしまいます．
5.実行結果
5.1.非正規化グラフラプラシアン
GClus = graph_clustering(K=15, n_clusters=5, lap_kind=None)の場合です．
確かに一つのノードが一つのクラスタになっていて，全然クラスタ感がしない感じになっちゃってますね．
5.2.酔歩正規化グラフラプラシアン
GClus = graph_clustering(K=15, n_clusters=5, lap_kind="rw")の場合です．
多少クラスタっぽくなってます．
5.3.クラスターの内訳
labels_split[0]
['MmeMagloire', 'CountessDeLo', 'Geborand']
labels_split[1]
['Myriel',
'Napoleon',
'MlleBaptistine',
'Champtercier',
'Cravatte',
'Count',
'OldMan',
'Labarre',
'Valjean',
'MmeDeR',
'Isabeau',
'Gervais',
'Tholomyes',
'Blacheville',
'Favourite',
'Dahlia',
'Zephine',
'MmeThenardier',
'Cosette',
'Javert',
'Fauchelevent',
'Bamatabois',
'Perpetue',
'Woman1',
'Judge',
'Champmathieu',
'Brevet',
'Chenildieu',
'Cochepaille',
'Pontmercy',
'Boulatruelle',
'Eponine',
'Anzelma',
'Woman2',
'MotherInnocent',
'Gribier',
'Jondrette',
'Gavroche',
'Gillenormand',
'Magnon',
'MlleGillenormand',
'MmePontmercy',
'MlleVaubois',
'LtGillenormand',
'Marius',
'BaronessT',
'Mabeuf',
'Enjolras',
'Courfeyrac',
'Bahorel',
'Bossuet',
'MotherPlutarch',
'Gueulemer',
'Babet',
'Montparnasse',
'Toussaint',
'Child2',
'MmeHucheloup']
labels_split[2]
['Combeferre', 'Prouvaire', 'Joly', 'Grantaire', 'Claquesous', 'Brujon']
labels_split[3]
['Marguerite',
'Listolier',
'Fameuil',
'Fantine',
'Thenardier',
'Simplice',
'MmeBurgon',
'Feuilly',
'Child1']
labels_split[4]
['Scaufflaire']
どのようなクラスタに別れたかというと，こんな感じです．僕はレミゼを映画でしか知らない上にジャンバルジャンしか覚えていないので，なんだかよくわからないです．ジャンバルジャンは一番大きいクラスタにいますね，しか言えない．