ラビットチャレンジ #10 サポートベクターマシン（SVM）

# サポートベクターマシン（SVM）
* 2クラス分類問題
    * 2クラス分類問題とは与えられた入力データが2つのカテゴリーのどちらに属する識別する問題

* SVMは2クラス分類問題の代表的な手法の一つで高い予測精度を持つ関数が構築できることが知られている
* SVMはもともと2クラス分類問題のために考案されたが、その後、回帰問題などへも応用されている

* このレポートにはSVMの2クラス分類問題について記述する。

### 2クラス分類
* 入力データを$x$という数値ベクトルで表す
* カテゴリAとカテゴリBに属するかを表す変数を$y$とする
    * 入力がAに分類されるときには$y=1$
    * 入力がBに分類されるときには$y=-1$をとるとする

* 2クラス問題とは何らかの入力（特徴量）$x$が与えられたときにラベル$y$の値を予想する問題
* つまり **$x$を与えたときに$y$を返す関数**を構築する！

### 決定関数と分類境界
* 入力$x$がどちらのクラスに属するか識別するために次の **決定関数**を使う

```math
f(x)=w^Tx+b
```
* 上記の式において$w$と$b$は事前にはわからない未知のパラメータであり、これを推定する
* このパラメータの推定方法は多数考案されており、 **SVMはその代表的な方法の1つ**

#### どのように分類されるか
* ある入力データ$x$に対して決定関数$f(x)$を計算し、その符号により2つのクラスに分類する

```math
y= sgn f(x) = \left\{
\begin{array}{ll}
+1 & (f(x) \gt 0) \\
-1 & (f(x) \lt 0)
\end{array}
\right.
```

* 変数$y$はラベル
* $sgn(・)$は符号関数
    * 引数が正の場合には+1を返す
    * 引数が負の場合には-1を返す
* つまり、決定関数$f(x)$の正負を調べることで入力$x$のクラスを識別できる

##### 特徴ベクトルが2次元のときを考えてみる
* 決定関数は下記のようになる

```math
f(x)=(w_1,w_2)\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} +b = w_1x_1+w_2x_2+b
```

* 上記の式は下記で表すような方程式となっている
    * $f(x) > 0$の場合（青の領域）では$y=1$
    * $f(x) < 0$の場合（赤の領域）では$y=-1$

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/4136e51b-e424-a1c6-728e-c5c14c923576.png)

* $x_1-x_2$平面と$f(x)$の平面の交戦である$f(x)=0$が特徴ベクトル$x={x_1, x_2}^T$をわける境界線となっているこれを **分類境界**という
* 式（決定関数）でパラメータ$w$と$b$を決定することは適切な分類教科を決定することに対応している

## 線形サポートベクター分類（ **ハードマージン**）
* SVMの種類の1つ
* n個の訓練データをすべて正しく分類できるパラメータ $w$と$b$がある場合
    * これを式で表すと $y_i f(x_i) > 0 (i=1,2,3,....,n)$となる
    * この式は訓練データの値（符号）$y_i$と予測値 $sign f(x_i)$の値（符号）が同じであることを表す
    * 正しく分類できていない場合は$y_i f(x_i) < 0$となる
    * 正しく分類できている式が成り立つ場合に、分類境界f(x) = 0を決定する方法を考えていく
    * 分類境界はいくつも考えられるがその中で最適なものを選ぶ
    * 分類境界から各クラスまでの距離を **マージン**という（下の図ではマージンは点線になる）
        * もう少し詳しく書くと赤クラスの分類境界に一番近いデータから青クラスの分類境界に1番近いデータまでの距離
    * 下のAとBは分類境界を引けているがマージンが異なる（Aのほうがマージンが小さく、Bのほうがマージンが大きい）
    * SVMではマージンが大きいほど良い分類境界となるため、この場合はBのほうが適切な分類境界
    * つまり、最適な分類境界はマージンが最大値をとる分類境界ということになる
    * この考え方を **マージン最大化**という
        * 分類境界 $f(x) = 0$と分類境界から最も近くにあるデータ$x_i$との距離の最大化
    * マージン最大化から、分類境界の決定に関わっているのは分類境界に最も近いデータ$x_i$のみで、他のデータは関係ないことがわかる
    * 実際、分類境界に最も近いデータ以外を取り除いてしまってもSVMの分類によって得られる分類境界は変化しないことが知られてる
    * 分類境界に最も近いデータ $x_i$が分類境界を支えていると解釈できるため、この $x_i$のことを **サポートベクトル**という

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/f2441122-6a5c-b9ac-5e54-ba5261cb76b2.png)

* この図のように赤と青が完全に分類できると仮定しているSVMのことを **ハードマージン**と呼ぶ

## 線形サポートベクター分類（ **ソフトマージン**）
* SVMの種類の1つ
* ハードマージンは訓練データを完璧に分類できる仮定していたが、訓練データが完璧に分類できないことを許容するのが **ソフトマージン**
* ソフトマージンは完璧に分類できるという仮定をなくし、SVMの分類を分離可能でないデータに適用できるように拡張していく
* ハードマージンからソフトマージンへの拡張は、制約条件を緩和し、多少の分類誤りは許すようにすることで実現する
* このために新たに非負の変数 $\varepsilon_i \geq 0  (i= 1,2,...,n)$を導入し16、上の制約条件を次のように変更

```math
y_i[w^Tx_i+b] \geq 1 - \varepsilon_i (i=1,....,n)
```

* この$\varepsilon_i$マージン内に入るデータや誤分類されたデータに対する誤差を表す変数でスラック変数と呼ばれる
* ハードマージンでは。マージンを分類境界とその分類境界に最も近い$x_i$との距離を用いて定義していたが、ソフトマージンでは下記図のように$f(x)=1$と$f(x)=-1$の間の距離をマージンという

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/3829cad4-3755-e2ef-24dd-39179d03fd03.png)

* ハードマージン同様、このマージンを最大化しつつ、分類の誤差$\varepsilon_i$を最小化するように分類境界を決定する


# 実装演習
## 線形分離可能な訓練データを作成する

```python
def gen_data():
    x0 = np.random.normal(size=50).reshape(-1, 2) - 2.
    x1 = np.random.normal(size=50).reshape(-1, 2) + 2.
    X_train = np.concatenate([x0, x1])
    ys_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
    return X_train, ys_train

X_train, ys_train = gen_data()
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/b1423108-6d8d-53f5-27ec-c4b2f372ad1d.png)

### 学習

```python
t = np.where(ys_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# 線形カーネル
K = X_train.dot(X_train.T)

eta1 = 0.01
eta2 = 0.001
n_iter = 500

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.where(a > 0, a, 0)
```

```python
index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()

xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)

# 訓練データを可視化
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域を可視化
#plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])


# マージンと決定境界を可視化
plt.quiver(0, 0, 0.1, 0.35, width=0.01, scale=1, color='pink')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/c4c028d3-788a-0792-2978-401c982798b9.png)

## 線形分離できない訓練データの作成

```python
factor = .2
n_samples = 50
linspace = np.linspace(0, 2 * np.pi, n_samples // 2 + 1)[:-1]
outer_circ_x = np.cos(linspace)
outer_circ_y = np.sin(linspace)
inner_circ_x = outer_circ_x * factor
inner_circ_y = outer_circ_y * factor

X = np.vstack((np.append(outer_circ_x, inner_circ_x),
               np.append(outer_circ_y, inner_circ_y))).T
y = np.hstack([np.zeros(n_samples // 2, dtype=np.intp),
               np.ones(n_samples // 2, dtype=np.intp)])
X += np.random.normal(scale=0.15, size=X.shape)
x_train = X
y_train = y

plt.scatter(x_train[:,0], x_train[:,1], c=y_train)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/f4285e26-9fce-faef-e6a6-f25c1bddd156.png)


### 学習

```python
def rbf(u, v):
        sigma = 0.8
        return np.exp(-0.5 * ((u - v)**2).sum() / sigma**2)
    
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# RBFカーネル
K = np.zeros((n_samples, n_samples))
for i in range(n_samples):
    for j in range(n_samples):
        K[i, j] = rbf(X_train[i], X_train[j])

eta1 = 0.01
eta2 = 0.001
n_iter = 5000

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.where(a > 0, a, 0)
```

### 予測

```python
index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()

xx0, xx1 = np.meshgrid(np.linspace(-1.5, 1.5, 100), np.linspace(-1.5, 1.5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * rbf(X_test[i], sv)
y_pred = np.sign(y_project)

# 訓練データを可視化
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域を可視化
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/45a31298-c308-76e7-a2e3-f1a5a1117b24.png)


## ソフトマージンSVM
* 完全に分離できないような訓練データを作成

```python
x0 = np.random.normal(size=50).reshape(-1, 2) - 1.
x1 = np.random.normal(size=50).reshape(-1, 2) + 1.
x_train = np.concatenate([x0, x1])
y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)

plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/cc54eaab-08a9-2a75-3cfb-888b0a5e304d.png)

### 学習

```python
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# 線形カーネル
K = X_train.dot(X_train.T)

C = 1
eta1 = 0.01
eta2 = 0.001
n_iter = 1000

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.clip(a, 0, C)
```

### 予測

```python
index = a > 1e-8
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()

xx0, xx1 = np.meshgrid(np.linspace(-4, 4, 100), np.linspace(-4, 4, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)

# 訓練データを可視化
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域を可視化
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2183979/c935e940-b698-e84d-a3d1-b8210f82d347.png)