はじめに
pythonでSVM（サポートベクターマシン）を実装しました．
教科書として『はじめてのパターン認識』を使いました．
本記事の構成
- はじめに
- SVM
- マージン最大化
- 主問題
- ラグランジュ未定乗数法
- KKT条件
- 双対問題
- 解の求め方
- pythonでの実装
- 結果
- おわりに
- 追記情報
SVM
SVMとは，最大マージンを実現する2クラス線形識別学習法です．
教科書によると，他クラスの識別は $K-1$ 個の一対多SVMで実現することが多いようです．
マージン最大化
次のようなデータ集合を考えます． データの総数を $N$ とすると，$i = 1,...,N$となります．
\begin{align}
& D_L = \bigl\{(t_i, {\boldsymbol{x}}_i)\bigr\} \\
& t_i = \{-1, +1\}
\end{align}
線形識別関数のマージンを $\kappa$ とすると，全ての学習データで次式が成り立ちます．
|\boldsymbol{w}^T{\boldsymbol{x}}_i + b| \geq \kappa
$\kappa$ で正規化したものを $\boldsymbol{w}$， $b$ と置き直し，教師との積をとることで絶対値をはずします．
t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) \geq 1 \tag{1}
図のようにクラス間マージンは，識別超平面に最も近いデータを $\boldsymbol{w}$ 方向に射影した長さの差として求めることができます．
$t_{i} = +1$ のクラスの中で識別超平面に最も近いデータは，$\boldsymbol{w}$ 方向に射影した長さがそのクラス内で最小となり，
$t_{i} = -1$ のクラスの中で識別超平面に最も近いデータは，$\boldsymbol{w}$ 方向に射影した長さがそのクラス内で最大となります．
クラス間マージン $\rho(\boldsymbol{w}, b)$ は次式で表されます．
\begin{align}
\rho(\boldsymbol{w}, b) &= \min_{\boldsymbol{x} \in C_{y=+1}}\frac{{\boldsymbol{w}}^T\boldsymbol{x}}{||\boldsymbol{w}||} - \max_{\boldsymbol{x} \in C_{y=-1}}\frac{{\boldsymbol{w}}^T\boldsymbol{x}}{||\boldsymbol{w}||} \\
\\
&= \frac{1-b}{||\boldsymbol{w}||} - \frac{-1-b}{||\boldsymbol{w}||} \\
\\
&= \frac{2}{||\boldsymbol{w}||} \tag{2}
\end{align}
式(2)より，$||\boldsymbol{w}||$ を最小化すれば，クラス間マージンは最大化されることが分かります．
また，クラス間マージンは識別超平面に最も近いデータだけで決まります．この識別超平面に最も近いデータのことをサポートベクターと呼びます．
主問題
式(1)の制約条件のもと，式(2)に表されるクラス間マージンを最大化することを考えます．これを主問題とします．
評価関数（最小化）
L_p(\boldsymbol{w}) = \frac{1}{2}{\boldsymbol{w}}^{T}\boldsymbol{w}
不等式制約条件
t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1 \geq 0
ラグランジュ未定乗数法
主問題を解くためにラグランジュ未定乗数法を使います．
ラグランジュの未定乗数法（ラグランジュのみていじょうすうほう、英: method of Lagrange multiplier）とは、束縛条件のもとで最適化を行うための数学（解析学）的な方法である。いくつかの変数に対して、いくつかの関数の値を固定するという束縛条件のもとで、別のある1つの関数の極値を求めるという問題を考える。各束縛条件に対して定数（未定乗数、Lagrange multiplier）を用意し、これらを係数とする線形結合を新しい関数（未定乗数も新たな変数とする）として考えることで、束縛問題を普通の極値問題として解くことができる方法である。
wikipediaより引用
具体例
何を言っているかさっぱり分からないので，具体例を挙げて説明します．
$g(x, y) = x^2 + y^2 - 2 = 0$ の条件下で $f(x, y) = y - x$ の最大化を考えます．
これは半径 $\sqrt{2}$ の円周上の $(x, y)$ を使って，直線の切片を最大化するという問題になります．
図から，円と直線が接するとき $f(x, y)$ が制約条件のもと最大値 $2$ を取ることが分かります．
これをラグランジュ未定乗数法を使って解いてみましょう．
\begin{align}
L(x, y, \lambda) &= f(x, y) + \lambda g(x, y) \\
\\
&= y - x + \lambda(x^2 + y^2 - 2) \tag{3}
\end{align}
式(3)のラグランジュ関数を $x, y, \lambda$ で偏微分し，それぞれイコール $0$ とします．
\begin{align}
& \frac{\partial L}{\partial x} = -1 + 2\lambda x = 0 \\
& \frac{\partial L}{\partial y} = 1 + 2\lambda y = 0 \\
& \frac{\partial L}{\partial \lambda} = x^2 + y^2 - 2 = 0
\end{align}
1，2番目の式を3番目の式に代入します．
\begin{align}
& \Bigl(\frac{1}{2\lambda}\Bigr)^2 + \Bigl(-\frac{1}{2\lambda}\Bigr)^2 - 2 = 0 \\
\\
& \lambda = \pm \frac{1}{2}, x = \pm 1, y = \mp 1
\end{align}
よって，$g(x, y) = 0$ の制約条件下で $f(x, y) = y - x$ は最大値 $2$ を取ります．
この最大値について3次元で考えてみます．3次元空間では，$z = f(x, y) = y - x$ は斜面に，$g(x, y) = x^2 + y^2 - 1 = 0$ は円柱になります．
$g(x, y) = 0$ の制約条件下で $z = y - x$ の最大値を求めたいので，斜面と円柱が交わる面に関して，$z$ 座標の最大値を見つければ良いこととなります．
つまり，$x^2 + y^2 - 2 = 0$ の円を $z = y - x$ 平面に射影したときの $z$ 軸方向の最大値が求めたい値です．
これを図で表します．真横から見たとき，斜面と円柱が交わっている部分の中で，$z$ 軸方向の最大値は $2$ であることが分かり，これは計算結果と一致します．
解説
なぜラグランジュ未定乗数法で，最小値を求められるのか解説してみます．（誤りを含むかもしれません）
制約条件 $g(\boldsymbol{x}) = 0$（以下，制約面とする）のもと，評価関数 $f(\boldsymbol{x})$ が $\boldsymbol{x}^\star$ において最大値をとるとき，$\boldsymbol{\nabla}f(\boldsymbol{x})$ は制約面に対して垂直でなければなりません．
なぜなら，垂直でないと仮定すると，制約面の方向に沿って $f(\boldsymbol{x})$ をさらに大きくすることができるからです．また，$\boldsymbol{\nabla}g(\boldsymbol{x})$ も制約面に対して垂直になります．
よって $\boldsymbol{x}^\star$ において $\boldsymbol{\nabla}f, \boldsymbol{\nabla}g$ は平行なベクトルとなり，次式が成り立ちます．
\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0 \tag{4}
ここで，ラグランジュ関数を以下のように定義します．
L(\boldsymbol{x}, \lambda) \equiv f(\boldsymbol{x}) + \lambda g(\boldsymbol{x})
ラグランジュ関数 $L(\boldsymbol{x}, \lambda)$ を $\boldsymbol{x}, \lambda$ で偏微分することで式(4)および制約条件を導くことができます．
これを解くことで停留点 $\boldsymbol{x}^\star$ が得られます．この $\boldsymbol{x}^\star$ の中に，$f(\boldsymbol{x})$ を最大化する $\boldsymbol{x}$ が存在します．
KKT条件
では，このラグランジュの未定乗数法を主問題に適用します．
ラグランジュ関数を次式のように定義します．
\tilde{L}_{p}(\boldsymbol{w}, b, \boldsymbol{\alpha}) = \frac{1}{2}{\boldsymbol{w}}^T\boldsymbol{w} - \sum_{i=1}^{N}\alpha_i(t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1) \tag{5}
ここで，$\boldsymbol{\alpha} = {(\alpha_1, ..., \alpha_N)}^T (\alpha_i \geq 0)$ であり，$\alpha_i$ がラグランジュ未定乗数です．
この最適化問題の解 $\boldsymbol{w}_0, b_0$ は，以下のKKT条件を満たす解として得られます．
< KKT条件 >
\begin{align}
& \left.\frac{\partial \tilde{L}_p(\boldsymbol{w}, b, \boldsymbol{\alpha})}{\partial \boldsymbol{w}}\right|_{\boldsymbol{w} = \boldsymbol{w}_0} = \boldsymbol{w}_0 - \sum_{i=1}^{N}\alpha_{i}t_i{\boldsymbol{x}}_i = 0 \tag{6} \\
& \frac{\partial \tilde{L}_p(\boldsymbol{w}, b, \boldsymbol{\alpha})}{\partial b} = \sum_{i=1}^{N}\alpha_{i}t_i = 0 \tag{7} \\
& \\
& t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1 \geq 0 \tag{8} \\
& \\
& \alpha_i \geq 0 \tag{9} \\
& \\
& \alpha_i(t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1) = 0 \tag{10}
\end{align}
式(6),(7)はラグランジュ関数の偏微分イコール $0$，式(8)は主問題の制約条件，式(10)は相補性条件です．
式(10)の相補性条件より
t_i(\boldsymbol{w}^T\boldsymbol{x}_i + b) - 1 > 0 \Rightarrow \alpha_i = 0 \\
t_i(\boldsymbol{w}^T\boldsymbol{x}_i + b) - 1 = 0 \Rightarrow \alpha_i \geq 0
相補性条件は，サポートベクター以外のデータによる制約条件の効果はなく、サポートベクターだけで識別超平面が決定されることを示しています．
（調べてみるとKKT条件はかなり奥が深そうなので，これくらいのふんわりとした理解で次に進むことにします）
双対問題
式(5)を展開し，最適解および式(6),(7)を代入します．
\begin{align}
L_{d}(\boldsymbol{\alpha}) &= \frac{1}{2}{\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} - \sum_{i=1}^{N}\alpha_{i}t_{i}{\boldsymbol{w}_{0}}^{T}{\boldsymbol{x}}_{i} - b\sum_{i=1}^{N}\alpha_{i}t_{i} + \sum_{i=1}^{N}\alpha_{i} \\
&= \frac{1}{2}{\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} - {\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} - 0 + \sum_{i=1}^{N}\alpha_{i} \\
&= \sum_{i=1}^{N}\alpha_{i} - \frac{1}{2}{\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} \\
&= \sum_{i=1}^{N}\alpha_{i} - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}t_{i}t_{j}{\boldsymbol{x}_i}^T{\boldsymbol{x}}_j
\end{align}
式(6)より，最適解がラグランジュ未定乗数と学習データの線形和で表されることから，$\alpha_{i}$ を求める問題に置き換えることができます．
最適なラグランジュ未定乗数は，$L_{d}(\boldsymbol{\alpha})$ を最大にする $\boldsymbol{\alpha}$ により得られます．これを主問題に対する双対問題と呼びます．
$N$ 個の $1$ を並べたベクトルを $\boldsymbol{1}$，学習データで作られた行列を $\boldsymbol{H}$，教師データベクトルを $\boldsymbol{t}$ とします．
\begin{align}
& \boldsymbol{1} = (1,...,1)^T \\
\\
& \boldsymbol{H} = (H_{ij} = {t_it_j\boldsymbol{x}_i}^T\boldsymbol{x}_j) \\
\\
& \boldsymbol{t} = (t_1,...,t_N)^T
\end{align}
評価関数（最大化）
L_{d}(\boldsymbol{\alpha}) = \boldsymbol{\alpha}^T\boldsymbol{1} - \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{H}\boldsymbol{\alpha}
等式制約条件
\boldsymbol{\alpha}^T\boldsymbol{t} = 0
双対問題を解いて得られたラングランジュ未定乗数の最適解を $\tilde{\boldsymbol{\alpha}} = (\tilde{\alpha}_1,...,\tilde{\alpha}_N)^T$ とすると，最適解は次のように求められます．
\begin{align}
& \boldsymbol{w}_0 = \sum_{i=1}^{N}\tilde{\alpha}_{i}t_i{\boldsymbol{x}}_i \\
& b_0 = \frac{1}{t_s} - {\boldsymbol{w}_0}^T\boldsymbol{x}_s
\end{align}
$\boldsymbol{x}_s, t_s$ はサポートベクターのデータおよび教師を表します．
解の求め方
主問題を双対問題に置き換え，最適な $\tilde{\boldsymbol{\alpha}}$ を求めれば，マージン最大化は実現できることを導きました．
$L_{d}(\boldsymbol{\alpha})$ を最大にする $\boldsymbol{\alpha}$ をプログラムで求めるために最急降下法を使います．
双対問題の制約条件を満たすよう罰金項を加え，以下のように $L_{d}(\boldsymbol{\alpha})$ を再定義しました．
L_{d}(\boldsymbol{\alpha}) = \boldsymbol{\alpha}^T\boldsymbol{1} - \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{H}\boldsymbol{\alpha} - \frac{1}{2}\beta\boldsymbol{\alpha}^T\boldsymbol{t}\boldsymbol{t}^T\boldsymbol{\alpha}
では，実際に最急降下法で $\alpha_i$ を更新します．
\begin{align}
{\alpha_i}' &= {\alpha_i} + \eta\frac{L_{d}(\boldsymbol{\alpha})}{\alpha_i} \\
&= {\alpha_i} + \eta(1 - \sum_{j=1}^{N}\alpha_{j}t_{i}t_{j}{\boldsymbol{x}_i}^T{\boldsymbol{x}}_j - \beta\sum_{j=1}^{N}\alpha_{j}t_{i}t_j)
\end{align}
ここまで難しい計算をしてきましたが，最終的に得られる数式は美しいですね．
これで $L_{d}(\boldsymbol{\alpha})$ を最大にする $\boldsymbol{\alpha}$ は得られます．
pythonでの実装
以下のように実装しました．
$\alpha_i$ の更新が簡単にできます．また，罰金項の効果で，$\boldsymbol{\alpha}^T\boldsymbol{t} = 0$ の条件を満たすように更新が行われます．
import numpy
from matplotlib import pyplot
import sys
def f(x, y):
return x - y
if __name__ == '__main__':
param = sys.argv
numpy.random.seed()
N = 30
d = 2
X = numpy.random.randn(N, d)
T = numpy.array([1 if f(x, y) > 0 else - 1 for x, y in X])
alpha = numpy.zeros(N)
beta = 1.0
eta_al = 0.0001 # update ratio of alpha
eta_be = 0.1 # update ratio of beta
itr = 1000
for _itr in range(itr):
for i in range(N):
delta = 1 - (T[i] * X[i]).dot(alpha * T * X.T).sum() - beta * T[i] * alpha.dot(T)
alpha[i] += eta_al * delta
for i in range(N):
beta += eta_be * alpha.dot(T) ** 2 / 2
index = alpha > 0
w = (alpha * T).T.dot(X)
b = (T[index] - X[index].dot(w)).mean()
if '-d' in param or '-s' in param:
seq = numpy.arange(-3, 3, 0.02)
pyplot.figure(figsize = (6, 6))
pyplot.xlim(-3, 3)
pyplot.ylim(-3, 3)
pyplot.plot(seq, -(w[0] * seq + b) / w[1], 'k-')
pyplot.plot(X[T == 1,0], X[T == 1,1], 'ro')
pyplot.plot(X[T == -1,0], X[T == -1,1], 'bo')
if '-s' in param:
pyplot.savefig('graph.png')
if '-d' in param:
pyplot.show()
結果
下の図のような結果が得られました．データ点の色がクラスを表しています．
計算により求められた直線を一緒に描画したところ，それっぽい結果は得られました．
おわりに
不等式条件下でマージン最大化を目的とした主問題を設定し，これをラグランジュ未定乗数法を使うことで双対問題へ置き換え，最適な $\boldsymbol{w}_0, b_0$ を求められることが分かりました．
追記情報
- 2016/05/01
@skyshk さんからコメントをいただいてコードを修正しました．