はじめに

以前、「機械学習の分類」で取り上げたアルゴリズムについて、その理論とpythonでの実装、scikit-learnを使った分析についてステップバイステップで学習していく。個人の学習用として書いてるので間違いなんかは大目に見て欲しいと思います。

今回はサポートベクターマシンについて。その歴史こそ古いものの、機械学習の分野では人気の手法です。汎化性能も高く、ILSVRC2012でディープラーニングに破られるまでは最強のアルゴリズムでした。これを覚えてるとだいぶよい気がします（語彙）。

今回参考にしたのは以下のサイト。ありがとうございます。他にも良書籍がたくさんあるようなのでそういうのもご参考に。

• PRML第7章のサポートベクターマシン(SVM)をPythonで実装
• SVM(サポートベクターマシン)
• サポートベクターマシンを理解する
• SVMについて簡単に理解する｡
• SVM(サポートベクターマシン)について
• サポートベクターマシンを手計算して理解する
• SMO徹底入門
• Python + NumPy で SMO 実装した

理論編

理論をしっかり理解しようと思いましたが、とても複雑で、おまけにこれ以上ないくらいしっかりと理論に触れているエントリがあったので、エッセンスだけ記述していこうと思います。

ざっくりと概論

サポートベクターマシン(Support Vector Machine)は、パーセプトロンやロジスティック回帰のような教師ありの二値分類器です。サポートベクターの考え方を導入することで、未知データに対する汎化性能が高く、カーネル法を用いることで線形分類が不可能な分類問題(XORなど)も分類できるようになります。
乱暴に言うとパーセプトロンにマージン最大化とカーネル法を持ち込んだ分類器と言えるようです。

理論編

次のような特徴量が2つに対し、二値分類が必要な状況を考えます。青い丸と赤い丸がそれぞれ学習データで、青い点を1(正例)、赤い点を-1(負例)とします。境界に引かれた緑の線を$$y=ax+b$$とします。正例、負例それぞれのうち、緑の線に一番近い点をサポートベクター、サポートベクターから境界までの距離をマージンと呼ぶことにします。

サポートベクターマシンは、教師データに対し、このマージンを最大化するような$a$と$b$を求めるという問題になります。
現実はこう都合よくいくケースばかりでもなく、きちんと分かれない場合が多くあるんですが、まずはこのシンプルなケースをもとに考えていきます。

サポートベクターマシンがいろんなケースでうまく分類できる理由は、このマージン最大化のおかげで、未知のデータに対してもいい感じに分類してくれるのと、境界面の近傍のデータを用いることで、外れ値にあまり影響しないという点があるようです。

初期設定

N個の教師データを$\boldsymbol{x}=(x_0, x_1, \cdots, x_{N-1})$、分類ラベルを $\boldsymbol{t}=(t_0, t_1, \cdots, t_{N-1})$、境界の式(識別関数という)を$g(x)=\boldsymbol{w}^Tx+w_0$とおきます。$\boldsymbol{w}$は、$\boldsymbol{x}$の重みベクトルです。

マージン最大化

$y=ax+b$とある点$x_n$との距離は$$\frac{|ax+b|}{\sqrt{a^2}}$$(参考)ですが、これを$g(x)$で考える。ある点$x_n$と$g(x)$の距離$|r|$は、$$|r|=\frac{|g(x)|}{|w|}$$と表せる。ラベル$t_n$は、$|t_n|=1$かつ、$t_ng(x)>=0$であることから、$$|r|=\frac{t_n(\boldsymbol{w}^Tx_n+w_0)}{|w|}$$となる。

$g(x)$から一番近い点$x_{min}$における$|r_{min}|$を最大化するために、$$|r_{min}|=\frac{t_{min}(\boldsymbol{w}^Tx_{min}+w_0)}{|w|}$$が最大となる$\boldsymbol{w}$と$w_0$を求めていきます。

制約条件

上記の式を最大化するには、分母が小さいほどよいですが、0になると不定になってしまう(解が一意に決まらない)ので、制約を加えます。
正例の点において、$g(x)=1$、負例は$g(x)=-1$であるとすると、$$t_ng(x)=t_n(\boldsymbol{w}^Tx_n+w_0) \geq 1$$という条件下で求めることにします。
これは、境界と最近傍の点$t_ng(x_{min})=1$とのマージン、$\frac{1}{|w|}$を最大化するので、$|w|$を最小化すればいいことになります。あとで微分することを考えて、$\frac{1}{2}|w|^2$を最小化すると置き換えても問題ありません。

整理すると、$$t_n(\boldsymbol{w}^Tx_n+w_0) \geq 1$$という条件下で、$$\frac{1}{2}|w|^2$$を最小化する$\boldsymbol{w}$と$w_0$を求めていきます。

ラグランジュの未定乗数法

上式のような制約条件のある関数を最小化するためには、ラグランジュの未定乗数法という手法を使います。これは、双対問題と言って、書き換えられた問題を解くことでもともとの問題を解いたことになるという理論を利用しています。

制約式$h(x)=t_n(\boldsymbol{w}^Tx_n+w_0)-1$のもとで、$f(x)=\frac{1}{2}|w|^2$を最大化するために、ラグランジュ乗数$\lambda$を導入したラグランジュ関数$$L(w, w_0, \lambda)=f(w, w_0)-\sum_{i=1}^{N}\lambda_ih(w, w_0)$$を定義します。

L(w,w_0,\lambda)=\frac{1}{2}|w|^2-\sum_{i=1}^{N}\lambda_i \{ t_i(\boldsymbol{w}^Tx_i+w_0)-1\}

なので、$w$と$w_0$について偏微分し、偏微分を0とおくとラグランジュ関数が$\lambda$だけの式になり、

L(\lambda)=\sum_{n=1}^{N}\lambda_n-\frac{1}{2}\sum_{n=1}^{N}\sum_{m=1}^{N}\lambda_n\lambda_mt_nt_mx_n^Tx_m

という式が導かれる。制約条件は、

\lambda_n\geq 0 \\
\sum_{i=1}^N\lambda_it_i=0

です。これで$L(\lambda)$を最大化する$\lambda$を求めればいいという問題に置き換えできたことになります。

上式を$w$で偏微分すると、
$$w=\sum_{i=1}^{N}\lambda_it_ix_i$$が求まります。

$w$が求まったので、$w_0$を求めます。$$g(x)=\boldsymbol{w}^Tx+w_0=\sum_{i=1}^{N}\lambda_it_ix_ix+w_0$$であり、サポートベクトルとのマージンを1とおいたので、全サポートベクトルについて、

t_n(\sum_{m} \lambda_mt_mx_mx_n+w_0)=1

が成り立ちます。実際には変形して

w_0=\frac{1}{N_M}\sum_{n}(t_n-\sum_{m}\lambda_mt_mx_mx_n)

として求めます。整理すると、$\lambda$が求めた後、

w=\sum_{i=1}^{N}\lambda_it_ix_i \\
w_0=\frac{1}{N_M}\sum_{n}(t_n-\sum_{m}\lambda_mt_mx_mx_n)

を最終的には計算で求めます。ではその$\lambda$はどのようにして求めればよいでしょうか。実際にはSMO(Sequential Minimal Optimization)という方法を使って求めることができます。

SMO

まず、解くべき問題を再掲します。

L(\lambda)=\sum_{n=1}^{N}\lambda_n-\frac{1}{2}\sum_{n=1}^{N}\sum_{m=1}^{N}\lambda_n\lambda_mt_nt_mx_n^Tx_m \\

制約条件

\lambda_n\geq 0 ,\sum_{i=1}^N\lambda_it_i=0

このように、制約条件のある$L(\lambda)$を最大化する問題は、$\lambda$をSMOを使って求めることができます。

Wikipediaによれば、SMOは日本語で逐次最小問題最適化法というそうです。勾配法のように反復しながら解に近づけていく手法で、任意の2変数を選択し、収束するまで繰り返します。
この任意の2変数のことをWorking Setと呼び、この2変数をどうやって選ぶかがキモになります。

流れとしては、

• KKT条件に反する変数$\lambda_1$を選択する
• $\lambda_2$を決める
• $\lambda_1, \lambda_2$を更新する

という処理をKKT条件に反する変数が存在しなくなるまで繰り返します。

KKT条件

カルーシュ・クーン・タッカー条件（Karush-Kuhn-Tucker condition）、以下KKT条件とは、一階導関数が満たすべき最適条件を指します。

実はKKT条件は上のラグランジュ未定乗数法のときにも使っており、

(1) $\frac{\partial{L}}{\partial{w}}=0$
(2) $\frac{\partial{L}}{\partial{w_0}}=0$
(3) $t_n(\boldsymbol{w}^Tx_n+w_0) \geq 1$
(4) $\lambda_n \geq 0$
(5) $\lambda_n(t_n(\boldsymbol{w}^Tx_n+w_0)-1) = 0$

という条件です。特に5番目の条件を相補性条件と言うそうです。

KKT条件違反のチェックして変数を1つ決定する

相補性条件$$\lambda_n(t_n(\boldsymbol{w}^Tx_n+w_0)-1) = 0$$を場合分けして、
(1) $\lambda_n=0$の場合、$t_n(\boldsymbol{w}^Tx_n+w_0) \geq 1$
(2) $\lambda_n > 0$の場合、 $t_n(\boldsymbol{w}^Tx_n+w_0) = 1$

が導かれるため、これを満たさない$\lambda$を選べばいいことになります。逆に$\lambda$が無くなったら解が求まったとします。ここで選んだ$\lambda$を$\lambda_1$とします。

変数をもう1つ決定する

もう一つの変数は、以下の順で選択します。まず、現在の$\lambda$で$\boldsymbol{w}$と$w_0$を求める必要があります。
この時の誤差関数を$$E_n=(\boldsymbol{w}^Tx_n+w_0)-t_n$$
とします。

(1) 変数の更新量が最大になるように選択

$\lambda_1$の場合の$E_1$との誤差が最大になるような$n$を選ぶ。

(2) 境界上に存在しない

$t_n(\boldsymbol{w}^Tx_n+w_0) = 1$ではない任意の点

(3) 残り

変数を更新する

更新する$\lambda_1$と$\lambda_2$は決定したが、更新する際には線形制約$$\sum_{i=1}^N\lambda_it_i=0$$があるため、この条件下で更新する必要がある。つまり、片方の$\lambda$を更新した場合は、もう一方の$\lambda$を調整しなければならない。更新後の$\lambda_1$と$\lambda_2$をそれぞれ$\lambda_1^{new}$、$\lambda_2^{new}$とすると

\lambda_1^{new}t_1+\lambda_2^{new}t_2=\lambda_1t_1+\lambda_2t_2

となります。これを$t_1=t_2$の場合と$t_1 \ne t_2$で場合分けし、

\lambda_1^{new}=\lambda_1+\frac{t_1(E_2-E_1)}{x_1^2+x_1x_2+x_2^2} \\
\lambda_2^{new}=\lambda_2+t_1t_2(\lambda_1-\lambda_1^{new})

を得ます。実際には、$\lambda_1$と$\lambda_2$の取りうる範囲を求めてクリッピングという処理が必要ですが力尽きました。

pythonの実装

まず最初に、理論編で力尽きたのと、長くなりすぎたのでscikit-learnの実装だけにします。ごめんなさい。いつか必ず自力で実装するつもりです。

scikit-learnの実装

pythonの実装はなんともあっさりです。scikit-learnでサポートベクトルマシンを使った分類をやるにはsklearn.svm.LinearSVCを使います。今回は分かりやすいように、正例50個負例50個がちゃんと分離されたサンプルを使います。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

%matplotlib inline

from sklearn import svm

fig, ax = plt.subplots()    

x1_1=np.ones(50)+10*np.random.random(50)
x1_2=np.ones(50)+10*np.random.random(50)
x2_1=-np.ones(50)-10*np.random.random(50)
x2_2=-np.ones(50)-10*np.random.random(50)

x1 = np.c_[x1_1,x1_2]
x2 = np.c_[x2_1,x2_2]
y  = np.array(np.r_[np.ones(50), -np.ones(50)])

model = svm.LinearSVC()
model.fit(np.array(np.r_[x1,x2]), y)

ax.scatter(x1[:,0],x1[:,1],marker='o',color='g',s=100)
ax.scatter(x2[:,0],x2[:,1],marker='s',color='b',s=100)

w = model.coef_[0]

x_fig = np.linspace(-12.,12.,100)
y_fig = [-w[0]/w[1]*xi-model.intercept_/w[1] for xi in x_fig]
ax.plot(x_fig, y_fig)

plt.show()

このように分類されました。あたり前ですね。

まとめ

最後力尽きてしまいましたが、雰囲気はつかめたんではないでしょうか（つかめてないですねw）。今回扱ったのは、基本中の基本、線形分離可能かつ正例と負例がきちんと分類できる例(ハードマージン)で導出しました。

より高度というか、現実に近いなサポートベクターマシンには、線形分離できないケース(カーネル法)と、正例と負例がきちんと分類できない例(ソフトマージン)を考える必要がありますが、基本的な理論は今回の内容とほとんど変わらないです。

次回以降ではその辺を扱っていきたいです。

追記:書きました