渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

パッケージユーザーのための機械学習(3)：サポートベクターマシン(SVM)

R Python 機械学習 サンプルデータで試す機械学習シリーズ

（※はてなフォトライフの不具合で正しくない順番で画像が表示されている可能性があります）


PythonでSMO-SVM書き下すという宿題がまだ終わってないくせにこれ書いていいのか物凄く迷うんですが（笑）、R Advent Calendar 2013の12月6日分第6回の担当に当たっているのでついでに書いちゃいます。


なのですが。実はその12月6日は米ネバダ州・タホ湖畔で開催中のNIPS 2013に参加中*1でupできるかどうか分からない*2ので、タイムスタンプ変えて予めupさせてもらいました。ルール破っちゃってごめんなさい。。。


ということで、今回の参考文献はこちら。未だに評は分かれるみたいですが、僕が推すのはいわゆるSVM赤本こと『サポートベクターマシン入門』です。


サポートベクターマシン入門

サポートベクターマシン入門


日本語訳がひどい*3とか色々言われていますが、それでもSVMに関するほぼ全てがこの本に書かれているので、特に原理的な特徴を知りたいのであればこれ以上のテキストはないと思います。


また、前にも書きましたがSMO (Sequential Maximal Optimization)による凸二次計画解法の実装法について触れた本はそもそも英文でも少なく、和書だとこれがほぼ唯一と言っても差し支えないので、特にアルゴリズム実装をやりたい人はPRMLに並んで持っておくべきだと個人的には思ってます。


そうそう、SMOまわりは最近実に素晴らしいslideshareが出たので是非こちらで共有を。


これは本当に大変貴重な資料です。実は赤本でもSMOの実装まわりの記述は不完全＆理論まわりも不親切なので、こういうしっかりまとまった資料があるとめちゃくちゃ助かります。


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


ひとまずSVMの雰囲気だけ見たいので、前回同様GitHubに置いてある以前のサンプルデータを使いましょう。コンバージョン(CV)に効くアクション(a1-a7)を探り出そうというテーマで用意されたデータです。これをdという名前でインポートしておきます。


ところで、RにはSVMを実行できるパッケージは複数あるんですが*4、今回は代表例として{e1071}のsvm()関数を使います。理由は簡単で、これはPython / Java / C++などの言語で使えるSVMのマルチプラットフォーム実装であるLIBSVMのRバージョンだからです。パラメータの与え方なども全く同じで、かなり細かくチューニングすることができます。


言い換えれば、これから見ていくsvm(){e1071}の特徴は全て他のPython / Java / C++などの言語でもLIBSVMを使うことでほぼ完全に再現可能だということです。こういうマルチプラットフォーム・パッケージはアドホック分析系にとってもアルゴリズム実装系にとっても非常に有用ですね。


ということで以下のような感じでやってみましょう。

> require("e1071")
 要求されたパッケージ e1071 をロード中です 
 要求されたパッケージ class をロード中です 
> d.svm<-svm(cv~.,d)
# formula式で与えるだけ、簡単

> summary(d.svm)

Call:
svm(formula = cv ~ ., data = d)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.1428571 

# デフォルト設定なのでガウシアンカーネル＆ソフトマージンSVM

Number of Support Vectors:  468

# サポートベクターは3000個中468個

 ( 216 252 )


Number of Classes:  2 

Levels: 
 No Yes

> table(d$cv,predict(d.svm,d[,-8]))
     
        No  Yes
  No  1402   98
  Yes   80 1420

# 元データに対する予測マトリクス。
# 正答率93.5%と、なかなか悪くない


典型的な機械学習分類器なので、こんな感じで素っ気なく結果が出ます（笑）。SVMは基本的にはそのままでは偏回帰係数とか変数重要度も出ないので*5、あくまでも分類結果での勝負になります。


SVMとは何ぞや


SVMと言えば誰もが知ってる機械学習の代表選手みたいな手法なので、こんなぐだぐだな記事で四の五の言うまでもないと思うんですが、一応書いておきます。SVMの本質は基本的には


の2点に尽きると思ってます。マージン最大化はそのままVC理論に支えられた汎化性能の高さ（≒分離超平面の滑らかさ）につながるし、カーネル法の導入はそれまでの線形識別関数では難しかった*6線形分離不可能パターンの正確な分離を可能にしています。即ち、「複雑に入り組んだデータであっても綺麗に分類できる上に未知データの分類性能も良い」機械学習分類器であると言えます。


ちょっと簡単にマージン最大化とカーネル法について触れておきましょう。マージン最大化自体は割と単純な発想で、下の図に見られる「マージン」\frac{b}{\|\mathbf{w}\|}を最大にするように、パラメータを決定するというものです。こうすることで、できるだけ2つのクラス間の距離が広く空くように分離超平面が定まります*7。


(wikipedia:en:Support_vector_machine)


この時、ラグランジュ未定乗数法で立式してこれを式変形すると凸二次最適化問題が同値の問題として得られるので、その最適解群\alpha_iをヒューリスティックに探せば良い、というのがSVMの基本です。その凸二次最適化問題は赤本にも載っているこいつです*8。

maximize W(\mathbf{\alpha})=\sum^{l}_{i=1} \alpha_i - \frac{1}{2} \sum^{l}_{i,j=1} y_i y_j \alpha_i \alpha_j K(\mathbf{x_i},\mathbf{x_j}),

subject to \sum^{l}_{i=1} y_i \alpha_i = 0, C \leq \alpha_i \leq 0, i = 1, \cdots, l.


ちなみにこれをバイアスを固定して2番目の式である線形制約をスキップし、単純なオンラインアルゴリズムに直した擬似コードが、赤本p.173に載っているこちら。

与えられたトレーニング集合Sと学習率\mathbf{\eta} \in (\mathbb{R}^{+})^lに対し

\mathbf{\alpha} \leftarrow \mathbf{0}

repeat

for i=1 to l

\alpha_i \leftarrow \alpha_i + \eta_i \( 1-y_i \sum^{l}_{j=1} \alpha_j y_j K(\mathbf{x}_i,\mathbf{x}j) \)

if \alpha_i < 0 then \alpha_i \leftarrow 0

else

if \alpha_i > C then \alpha_i \leftarrow C

end for

until 停止基準を満たす

return \mathbf{\alpha}


一応これを組んでも全ての\alpha_iが求まりSVMっぽい分類器が出来上がりますが、線形制約が満たされていないので微妙に汎化性能の悪い分離超平面を描く羽目になります*9。


なお、このアルゴリズムだと下のif文のところで非ゼロとなる\alpha_iを選別しているので、学習データ全体ではなく一部のデータ（サポートベクター）のみが学習に関与することになります。これがいわゆるSVMのスパース性ってやつですね。サポートベクターが少なければ少ないほど効率よく、かつ汎化性能が高く学習できて、オーバーフィッティングを抑制できます*10。


ここでバイアスも含めて最適解に到達するよう線形制約込みでヒューリスティックに凸二次計画を解くアルゴリズムがSMOで、これについては赤本を上回る素晴らしい資料が上記の通りslideshareに出ている（『SMO徹底入門 - SVMをちゃんと実装する』）ので是非そちらをご参照あれ。


一方、カーネル法は一言でいえば「適当な空間写像を作ることで元の空間上では線形分離不可能なものを無理やり線形分離可能に変えてしまう」というもので*11、高次元特徴空間における内積計算をシンプルな形に直せてしまうというその数学的性質から「カーネルトリック」とも呼ばれます。


ちなみにカーネル関数には色々な種類があり、線形内積カーネル、RBFカーネル（ガウシアンカーネル*12）、多項式カーネル、シグモイドカーネルなどを適宜選んで、性能比較することでより良いものを選ぶのが常道です。


ちなみに多項式カーネルについてはその様子をCGで可視化した動画がYouTubeにあって見ることができます。



元々線形分離不可能だったデータが、高次元空間上では線形分離可能パターンに置き換わっていて、スッパリと分離超平面で分類できることが分かりやすく見て取れると思います。


決定境界を描いてみる


それではやってみましょう。いきなりですが、SVMの威力を確かめるためにXOR（線形分離不可能）パターンでやってみようと思います。GitHubからXORパターンのシンプル版、複雑版を持ってきて、それぞれxors, xorcという名前でインポートしておきましょう。

> require("e1071")
 要求されたパッケージ e1071 をロード中です 
 要求されたパッケージ class をロード中です 

> xors$label<-as.factor(xors$label)
> xorc$label<-as.factor(xorc$label)
# 学習ラベルをfactor型に直す

> xors.svm<-svm(label~.,xors)
> xorc.svm<-svm(label~.,xorc)
# それぞれSVMで推定する
# デフォルトではガウシアンカーネルが選択される

> px<-seq(-3,3,0.03)
> py<-seq(-3,3,0.03)
> pgrid<-expand.grid(px,py)
> names(pgrid)<-c("x","y")
# 分離超平面描画用のグリッド

> out.xors.svm<-predict(xors.svm,pgrid,type="vector")
> out.xorc.svm<-predict(xorc.svm,pgrid,type="vector")
# グリッド値をpredictで与える

> plot(xors[1:50,-3],col="blue",pch=19,cex=3,xlim=c(-3,3),ylim=c(-3,3))
> points(xors[51:100,-3],col="red",pch=19,cex=3)
> par(new=T)
> contour(px,py,array(out.xors.svm,dim=c(length(px),length(py))),xlim=c(-3,3),ylim=c(-3,3),col="purple",lwd=3,drawlabels=F)
# シンプルパターン

> plot(xorc[1:50,-3],col="blue",pch=19,cex=3,xlim=c(-3,3),ylim=c(-3,3))
> points(xorc[51:100,-3],col="red",pch=19,cex=3)
> par(new=T)
> contour(px,py,array(out.xorc.svm,dim=c(length(px),length(py))),xlim=c(-3,3),ylim=c(-3,3),col="purple",lwd=3,drawlabels=F)
# 複雑パターン

> table(xors$label,predict(xors.svm,xors[,-3]))
   
     1  2
  1 48  2
  2  0 50

# シンプルパターンの予測マトリクス

> table(xorc$label,predict(xorc.svm,xorc[,-3]))
   
     1  2
  1 36 14
  2  6 44

# 複雑パターンの予測マトリクス


f:id:TJO:20131202165748p:plain

シンプルパターンの場合と、

f:id:TJO:20131229154923p:plain

複雑パターンの場合。


シンプルパターンでは98%の正答率とまずまず。SVMらしい汎化性能がよく出ていて、非常に滑らかで美しい分離超平面を描いています。これは決定木の時とは大きな違いです。


一方、複雑パターンでは80%しか正答率が出てません。汎化が強過ぎて、明らかに真ん中のあたりが分類エラーだらけになってます。これではいくら美しい分離超平面を描いていても、機械学習分類器としては良くないですね。


そこで、ガウシアンカーネルで使えるハイパーパラメータであるcost*13とgamma*14をいじって、極端に汎化を抑えてオーバーフィッティング気味にしてみるとこうなります。

> xorc.svm2<-svm(label~.,xorc,cost=50,gamma=10)
> out.xorc.svm2<-predict(xorc.svm2,pgrid,type="vector")
> plot(xorc[1:50,-3],col="blue",pch=19,cex=3,xlim=c(-3,3),ylim=c(-3,3))
> points(xorc[51:100,-3],col="red",pch=19,cex=3)
> par(new=T)
> contour(px,py,array(out.xorc.svm2,dim=c(length(px),length(py))),xlim=c(-3,3),ylim=c(-3,3),col="purple",lwd=3,drawlabels=F)
> table(xorc$label,predict(xorc.svm2,xorc[,-3]))
   
     1  2
  1 50  0
  2  1 49


f:id:TJO:20131202170212p:plain


正答率99%まで上がったけど、これじゃSVMやる意味ないじゃん（笑）。あんまりなので、もうちょっと緩めてみます。

> xorc.svm3<-svm(label~.,xorc,cost=15,gamma=5)
> out.xorc.svm3<-predict(xorc.svm3,pgrid,type="vector")
> plot(xorc[1:50,-3],col="blue",pch=19,cex=3,xlim=c(-3,3),ylim=c(-3,3))
> points(xorc[51:100,-3],col="red",pch=19,cex=3)
> par(new=T)
> contour(px,py,array(out.xorc.svm3,dim=c(length(px),length(py))),xlim=c(-3,3),ylim=c(-3,3),col="purple",lwd=3,drawlabels=F)
> table(xorc$label,predict(xorc.svm3,xorc[,-3]))
   
     1  2
  1 45  5
  2  4 46


f:id:TJO:20131202171027p:plain


正答率91%になりましたが、全然汎化できてる気がしないですね。。。ともあれ、ぐちゃぐちゃに入り組んだXORパターンでもSVMなら（やろうと思えば）100%に近い正答率を叩き出すことができる、ということはこれで分かりました。


注意点


カーネルパラメータの設定と、カーネルそのものの選択が最大のポイントです。まず、カーネルパラメータを変えることで分類性能が大きく変わるのが普通なので、グリッドサーチなどでできるだけパラメータの最適化を図るべきです。


実際問題、上の例でも見た通りcostとgammaをちょっと変えるだけでかなり分類性能が変わってきますし、そもそも分離超平面の形状そのものが大きく変わります。この辺はグリッドサーチをかけるだけでなく、データの「非線型性」を感じ取る「センス」も重要になってくるかもしれません。


そして、カーネルそのものの選択について。LIBSVMではデフォルトだとガウシアンカーネルが選択されますが、ガウシアンカーネルだからといって必ずしも万能なわけではありません。線形分離可能パターンに対してガウシアンカーネルを選ぶと、逆に分類性能が下がるケースが存在することが知られています。


また、単に計算負荷の点で線形分離可能でSVMではなくても十分イケるものに対して、わざわざSVMをかけることで余計なリソースを食う可能性もあります。例えばロジスティック回帰の時に用いた二値分類データに対して用いると（dlという名前でインポートしておく）、

> dl$label<-as.factor(dl$label)

> ps<-seq(0,4,0.03)
> pt<-seq(0,5,0.03)
> pgrid2<-expand.grid(ps,pt)
> names(pgrid2)<-c("x","y")

# ガウシアンカーネル非線形SVM
> out.dl.svm1<-predict(dl.svm1,pgrid2,type="vector")
> plot(dl[1:25,-3],col="blue",pch=19,cex=3,xlim=c(0,4),ylim=c(0,5))
> points(dl[26:50,-3],col="red",pch=19,cex=3)
> par(new=T)
> contour(ps,pt,array(out.dl.svm1,dim=c(length(ps),length(pt))),xlim=c(0,4),ylim=c(0,5),col="purple",lwd=3,drawlabels=F)
> table(dl$label,predict(dl.svm1,dl[,-3]))
   
     0  1
  0 22  3
  1  2 23

# 前回取り上げたロジスティック回帰
> out.glm<-predict(dl.glm,pgrid2)
> plot(dl[1:25,-3],pch=19,cex=3,col="red",xlim=c(0,4),ylim=c(0,5))
> points(dl[26:50,-3],pch=19,cex=3,col="blue")
> par(new=T)
> contour(ps,pt,array(out.glm,dim=c(length(ps),length(pt))),xlim=c(0,4),ylim=c(0,5),col="purple",lwd=3,drawlabels=T,levels=0.5)
> table(dl$label,(sign(predict(dl.glm,dl[,-3]))+1)/2)
   
     0  1
  0 22  3
  1  3 22


f:id:TJO:20131202173019p:plain

ガウシアンカーネル非線形SVMの場合と、

f:id:TJO:20131203120044p:plain

前回取り上げたロジスティック回帰の場合。


見ての通り、たった1%差ということでほとんど分類性能に差が出ません。


特にスパム判定など、素性がほとんど互いに重ならないもの同士で分類するようなケースでは多次元空間上では線形分離可能なことが多く、注意が必要です。


最後に


SVMはVC理論に支えられた汎化性能の高さゆえ分離超平面が美しいので、こうして描画してみると楽しいんですよね。それに引き換え、他の分類器と来たら。。。おっと（笑）。

*1:うちの教授に引っ張られて行きます（笑）

*2:ネットワークの関係などで

*3:誤訳が多いとか文章が乱れているとか

*4:{kernlab}や{caret}など

*5:モンテカルロ法を使えば出せますが、手間ですよね

*6:逆伝播法ではできますが

*7:ただし誤分類が生じるのを前提としていて、VC理論によってその上限値が与えられるようになっています

*8:ただし1-ノルムソルトマージンの場合

*9:やってみると分かりますが少しカクカクします

*10:赤本p.138にVC理論に基づく汎化性能の定理が示されています

*11:みんな大好き再生核ヒルベルト空間のお話ですよ！

*12:厳密にはRBFカーネルはガウシアンだけではないが、普通はガウシアンでOK

*13:1-ノルムソフトマージンのC

*14:ガウシアンカーネルの式がsvm(){e1071}ではexp(-gamma*|u-v|^2)と書かれているので、そのgamma。つまりカーネルの幅と思えば大体合ってる