pythonでSVM実装

# はじめに

pythonでSVM（サポートベクターマシン）を実装しました．
教科書として[『はじめてのパターン認識』](http://www.amazon.co.jp/はじめてのパターン認識-平井-有三/dp/4627849710)を使いました．

本記事の構成

* はじめに
* SVM
    * マージン最大化
    * 主問題
    * ラグランジュ未定乗数法
    * KKT条件
    * 双対問題
* 解の求め方
* pythonでの実装
* 結果
* おわりに
* 追記情報


# SVM

SVMとは，最大マージンを実現する2クラス線形識別学習法です．
教科書によると，他クラスの識別は $K-1$ 個の一対多SVMで実現することが多いようです．

## マージン最大化

次のようなデータ集合を考えます． データの総数を $N$ とすると，$i = 1,...,N$となります．

```math
\begin{align}
& D_L = \bigl\{(t_i, {\boldsymbol{x}}_i)\bigr\} \\
& t_i = \{-1, +1\}
\end{align}
```

線形識別関数のマージンを $\kappa$ とすると，全ての学習データで次式が成り立ちます．

```math
|\boldsymbol{w}^T{\boldsymbol{x}}_i + b| \geq \kappa
```

$\kappa$ で正規化したものを $\boldsymbol{w}$， $b$ と置き直し，教師との積をとることで絶対値をはずします．

```math
t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) \geq 1 \tag{1}
```

図のようにクラス間マージンは，識別超平面に最も近いデータを $\boldsymbol{w}$ 方向に射影した長さの差として求めることができます．

<img width=420 alt="mergin.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/7bb67956-8b00-fa8a-c832-52fd8072aff6.png">

$t_{i} = +1$ のクラスの中で識別超平面に最も近いデータは，$\boldsymbol{w}$ 方向に射影した長さがそのクラス内で最小となり，
$t_{i} = -1$ のクラスの中で識別超平面に最も近いデータは，$\boldsymbol{w}$ 方向に射影した長さがそのクラス内で最大となります．
クラス間マージン $\rho(\boldsymbol{w}, b)$ は次式で表されます．

```math
\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)に表されるクラス間マージンを最大化することを考えます．これを主問題とします．

評価関数（最小化）

```math
L_p(\boldsymbol{w}) = \frac{1}{2}{\boldsymbol{w}}^{T}\boldsymbol{w}
```

不等式制約条件

```math
t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1 \geq 0
```

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

主問題を解くためにラグランジュ未定乗数法を使います．
> ラグランジュの未定乗数法（ラグランジュのみていじょうすうほう、英: method of Lagrange multiplier）とは、束縛条件のもとで最適化を行うための数学（解析学）的な方法である。いくつかの変数に対して、いくつかの関数の値を固定するという束縛条件のもとで、別のある1つの関数の極値を求めるという問題を考える。各束縛条件に対して定数（未定乗数、Lagrange multiplier）を用意し、これらを係数とする線形結合を新しい関数（未定乗数も新たな変数とする）として考えることで、束縛問題を普通の極値問題として解くことができる方法である。
[wikipedia](https://ja.wikipedia.org/wiki/ラグランジュの未定乗数法)より引用

### 具体例
何を言っているかさっぱり分からないので，具体例を挙げて説明します．
$g(x, y) = x^2 + y^2 - 2 = 0$ の条件下で $f(x, y) = y - x$ の最大化を考えます．
これは半径 $\sqrt{2}$ の円周上の $(x, y)$ を使って，直線の切片を最大化するという問題になります．

<img width=360 alt="rag2.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/a4c5207b-7bb5-c525-d5bb-1a851abdc91d.png">

図から，円と直線が接するとき $f(x, y)$ が制約条件のもと最大値 $2$ を取ることが分かります．
これをラグランジュ未定乗数法を使って解いてみましょう．

```math
\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$ とします．

```math
\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番目の式に代入します．

```math
\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$ であることが分かり，これは計算結果と一致します．

<img width=600 alt="lagrange.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/aca4e678-adb7-0689-2da9-85e032ea43e0.png">

### 解説
なぜラグランジュ未定乗数法で，最小値を求められるのか解説してみます．（誤りを含むかもしれません）
制約条件 $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$ は平行なベクトルとなり，次式が成り立ちます．

```math
\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0 \tag{4}
```

ここで，ラグランジュ関数を以下のように定義します．

```math
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条件

では，このラグランジュの未定乗数法を主問題に適用します．
ラグランジュ関数を次式のように定義します．

```math
\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条件 >**

```math
\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)の相補性条件より

```math
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)を代入します．

```math
\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}$ とします．

```math
\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}
```

評価関数（最大化）

```math
L_{d}(\boldsymbol{\alpha}) = \boldsymbol{\alpha}^T\boldsymbol{1} - \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{H}\boldsymbol{\alpha}
```

等式制約条件

```math
\boldsymbol{\alpha}^T\boldsymbol{t} = 0
```

双対問題を解いて得られたラングランジュ未定乗数の最適解を $\tilde{\boldsymbol{\alpha}} = (\tilde{\alpha}_1,...,\tilde{\alpha}_N)^T$ とすると，最適解は次のように求められます．

```math
\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})$ を再定義しました．

```math
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$ を更新します．

```math
\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$ の条件を満たすように更新が行われます．

```svm.py
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()
```


# 結果

下の図のような結果が得られました．データ点の色がクラスを表しています．
計算により求められた直線を一緒に描画したところ，それっぽい結果は得られました．

<img width=420 alt="result.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/d38d90f8-7a4f-e473-317a-652a69bf477c.png">


# おわりに

不等式条件下でマージン最大化を目的とした主問題を設定し，これをラグランジュ未定乗数法を使うことで双対問題へ置き換え，最適な $\boldsymbol{w}_0, b_0$ を求められることが分かりました．


# 追記情報

* 2016/05/01
@skyshk さんからコメントをいただいてコードを修正しました．