パッケージユーザーのための機械学習(10)：Affinity Propagation - 渋谷駅前で働くデータサイエンティストのブログ

だいぶ前回から間が空いてしまいましたが、ついに10回目になったこのシリーズ記事。。。多分クラスタリングだとこれが最後になるんじゃないでしょうか。以前話題に出ていたAffinity Propagationをやってみようと思います。

なのですが。今回も文献資料は見つけられずorz　ということでweb上の資料を適宜参照しながらやっていく形でいってみようかと思います。今回参考にさせていただいたのはこちら。

一番上はこのAffinity Propagationのアルゴリズムを考案した張本人であるFrey博士のサイト、次がググるとヒットするnojima718さんのブログ記事、最後がそのブログ記事中で紹介されているどこかの大学の方の輪講資料（勝手にご紹介させていただいてすみません）のPDFです。この3つを読めば、何となくどんなものかは分かるのではないかなぁと思います。

まずRでどんなものか見てみる

めっきり教師なし学習になってから微妙な感じのあるこのコーナーですが、めげずにやってみます。GitHubに置いてある多変量データを落としてきて、dという名前でインポートしておきます。Affinity PropagationはCRANパッケージでは{apcluster}で実装されています。

> library("apcluster", lib.loc="C:/Program Files/R/R-3.0.2/library")
 要求されたパッケージ Rcpp をロード中です 

 次のパッケージを付け加えます: ‘apcluster’ 

 以下のオブジェクトはマスクされています (from ‘package:stats’) : 

     heatmap 

> d.apc<-apclusterL(s=negDistMat(r=2),x=d[,-8],frac=0.02,sweeps=2,p=-7)
> d.apc.hcl<-as.hclust(d.apc)
> d.apc.hcl.cut<-cutree(d.apc.hcl,k=2)
> d.apc.labels<-labels(d.apc,type="enum")
> tmp<-rep(0,3000)
> for (i in 1:3000) {
+ tmp[i]<-d.apc.hcl.cut[d.apc.labels[i]]
+ }
> table(tmp,d$cv)
   
tmp   No  Yes
  1 1289  234
  2  211 1266

# 正答率85.1%。。。低いな。。。

思ったよりパフォーマンス良くないですねorz　ちょっと多次元＆（それなりの）大容量データでどうしたらいいかは僕も分かってないのでまた改めて。。。

Affinity Propagationとは何ぞや

ぶっちゃけ冒頭に挙げた資料を読んでいただければ十分です、ってかそれ以上の知識は僕にもありませんorz　なので資料から引用してちょっと僕の分かる範囲でかいつまんで以下にまとめてみました。

Affinity Propagation(AP) では全ての要素 (データ点)が exemplar になるかもしれないということを考慮し，exemplar とクラスタが良い組み合わせになるようにメッセージを交換する．そのため，各要素は常に「クラスタ中心になる」「クラスタ (の中心) に属する」という 2種類の面が考慮される．全要素間の関係性を similarity として設定し，それを元に exemplar 候補からクラスタメンバー (クラスタに属する点) に送られる responsibility とクラスタメンバーから exemplar 候補に送られる availability という2種類のメッセージ (一種の評価値) を交換しあう．各要素を i，exemplar 候補を k とおいて，responsibility とavailability をそれぞれ r(i, k)，a(i, k) と表現する．各要素間でメッセージの交換を繰り返して最終的に a(i, k) +r(i, k) で評価し，exemplar を決定，クラスタを生成する．
http://stanag4172.m13.coreserver.jp/document/AffinityPropagation.rev2.pdf

基本的には初期値依存性がないという点でK-meansよりも優れているというのがAffinity Propagationのセールスポイントで、初期値を変えて繰り返すという手間が要らない分計算時間がK-meansよりも短いということのようです。

上記引用の通り、APで重要なのはresponsibilityとavailabilityという2つの評価値。資料によれば、similarityを含めたそれぞれの定義は以下の通りということだそうです。

similarity
$s(i,k)$
負のユークリッド距離（負の二乗誤差）などが用いられる。Rでは他にも色々利用できる
responsibility
$r(i,k) \gets s(i,k) - \max_{k' s.t. k' \neq k} \{ a(i,k') + r(i,k') \}$
（ただし $k'$ は競合しているexemplar候補、メッセージの計算はresponsibilityから始めるが、その際に $a(i,k)=0$ と初期化される）
availability
$a(i,k) \gets \min \{ 0, r(k,i) + \displaystyle{\sum_{i' s.t. i' \not\in i,k} \max \{ 0, r(i',i) \} } \}$
（ただし $i'$ はexemplar候補 $k$ を最も高く評価しているデータ点）

これらのデータと2種類の評価値を用いて、資料中Fig.6のようなフローチャートに基づく繰り返し計算アルゴリズムによって学習していくというのが、APの基本形ということみたいですねー。

真の分布の分かっているデータをクラスタリングしてみる

そんなわけで、いつも通りRで試してみましょう。GitHubからシンプルパターン、複雑パターンそれぞれを落としてきてxors, xorcという名前でインポートしておきます。これらのデータの正体は以下の通りです。

# シンプルなXORパターン
> p11<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p12<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p13<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> p14<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> xors<-as.data.frame(cbind(rbind(p11,p13,p12,p14),t))
> names(xors)<-c("x","y","label")

# 複雑なXORパターン
> p21<-cbind(rnorm(n=25,mean=1,sd=1),rnorm(n=25,mean=1,sd=1))
> p22<-cbind(rnorm(n=25,mean=-1,sd=1),rnorm(n=25,mean=1,sd=1))
> p23<-cbind(rnorm(n=25,mean=-1,sd=1),rnorm(n=25,mean=-1,sd=1))
> p24<-cbind(rnorm(n=25,mean=1,sd=1),rnorm(n=25,mean=-1,sd=1))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> xorc<-as.data.frame(cbind(rbind(p21,p23,p22,p24),t))
> names(xorc)<-c("x","y","label")

(1,1), (-1,1), (-1,-1), (1,-1)のそれぞれを中心とする標準偏差0.5ないし1.0の二次元正規分布が4つ配置されたデータですね。ということで、xors / xorc双方のデータに対して{apcluster}パッケージを用いてAffinity Propagationによるクラスタリングを行ってみましょう。

ちなみに通常はapcluster()関数でやれますが、データ行数が多い場合はapclusterL()関数でデータを分割しながらやる必要が出てきます（冒頭の多変量データでの実行例を参照のこと）。クラスタリング結果はplot()メソッドを使うことで、exemplarと当該クラスタ所属サンプルとを明示的にプロットすることができます。その実行結果は以下の通り。

> require("apcluster")

> xors.apc<-apcluster(s=negDistMat(x=xors,r=2),q=0.1)
> plot(xors.apc,xors)

> xorc.apc<-apcluster(s=negDistMat(x=xorc,r=2),q=0.025)
> plot(xorc.apc,xorc)

f:id:TJO:20140728184522p:plain

シンプルパターンの結果と、

f:id:TJO:20140728184531p:plain

複雑パターンの結果。前者はともかく、後者は明らかにexemplarの位置が真の分布の中心からかなり外れているので、推定結果としては微妙かなぁ。。。

ちなみに何となく試行錯誤しながらチューニングして、この形に持ってきました。なので、自動ではこの結果にならないと思います（汗）。そういう意味でいうと、正直言ってパラメータチューニングをほとんどしなくても初期値さえちょこちょこいじれば概ねまともな結果になるK-meansに比べるとちょっと扱いづらいかなぁという気がします。

後はクラスタ数がバンバン変化する点ですね（やってみるとよく分かるはず）。これも混合ディリクレ過程みたいな感じで、どうやって使ったものかなぁみたいな部分があるので色々悩ましい点でもあり。。。ちょっとまた勉強してから考えてみます。。。