numpyとNetworkXで関係データ学習（スペクトラルクラスタリング）

# 1.はじめに
この記事は， **pythonによるネットワーク分析（コミュニティ抽出）**，ネットワーク分析に便利なライブラリである **NetworkXの簡単な使い方のまとめ**です．
コミュニティ抽出のアルゴリズムとして， **スペクトラルクラスタリングをnumpyで実装**します．


個人的なモチベーションとしては，この記事（[Probabilistic Matrix Factorization を導出して Edward で実装する](http://yamaguchiyuto.hatenablog.com/entry/2017/07/13/080000)）を読んで，この本（[関係データ学習](https://www.amazon.co.jp/dp/4061529218)）面白そうだなと思って読み始めたのでそのまとめのようなものです．

# 2.環境構築
## 2.1.今回の実施環境
- Python 3.6.0 :: Anaconda 4.3.1 (x86_64)
- numpy 1.11.3

## 2.2.インストール
```sh
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.データのソース
[ここ](http://www-personal.umich.edu/~mejn/netdata/)にネットワーク分析関係のお試しデータがいくつも公開されているので，この中から，


> 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)：『レ・ミゼラブル』の人物相関図](http://tjo.hatenablog.com/entry/2015/04/10/080000)

同じコミュニティ抽出の課題を，別の手法で試しています．このレミゼデータがどんなデータで，コミュニティ抽出をするとどんなことがわかるのかというのはこの記事を読むと大体わかるのではと思います．

この記事を読んだ感じ，ネットワーク分析のライブラリはpythonよりもRの方が揃っているのではと思いました．もしくは[cytoscape](http://www.cytoscape.org/)などを使うべきなのでしょうか．（もしくは僕が他のものを知らないだけ）

## 3.3.データの形式
```.gml```という拡張子のデータが，上記リンクのzipの中に入っています．gmlとは，Graph Modelling Languageの略で，グラフ構造のデータの保存形式です．

以下のような構造で，グラフのノードとエッジの情報を格納しています．

```sh
$ 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.クラスや関数の定義
実際にスペクトラルクラスタリングを行うクラスと，実行結果をプロットするのに必要な関数を定義します．

```python
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.コードを動かす
```python
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)```の場合です．
確かに一つのノードが一つのクラスタになっていて，全然クラスタ感がしない感じになっちゃってますね．

![fig2.png](https://qiita-image-store.s3.amazonaws.com/0/164264/cfc6bb94-03df-6692-5b01-97cfd09a3e37.png "fig2.png")




## 5.2.酔歩正規化グラフラプラシアン
```GClus = graph_clustering(K=15, n_clusters=5, lap_kind="rw")```の場合です．
多少クラスタっぽくなってます．

![fig1.png](https://qiita-image-store.s3.amazonaws.com/0/164264/fbecdc08-201f-a0d0-6784-d6f219e17855.png "fig1.png")

## 5.3.クラスターの内訳
```python
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']
```

どのようなクラスタに別れたかというと，こんな感じです．僕はレミゼを映画でしか知らない上にジャンバルジャンしか覚えていないので，なんだかよくわからないです．ジャンバルジャンは一番大きいクラスタにいますね，しか言えない．


# 6.終わりに
- 理論的な話については今の僕の知識では[関係データ学習](https://www.amazon.co.jp/dp/4061529218)をコピペするだけになってしまって面倒and意味なしなので省略して，この本を参照してください．いい本です．
- ちなみにこの本にはあともう一つ対称正規化ラプラシアンによるクラスタリングも登場します．クラスタリングはないですが，対称正規化ラプラシアンの計算はNetworkXに関数が定義されていて，[ここ](https://networkx.github.io/documentation/networkx-1.10/reference/generated/networkx.linalg.laplacianmatrix.normalized_laplacian_matrix.html#networkx.linalg.laplacianmatrix.normalized_laplacian_matrix)に詳細があります．
- クラスタリングを実施するときに，固有値の小さい方から着目するのですが，それが新鮮で面白いです．PCAとかだと大きい方に着目しますし．