その他

E資格 学習内容まとめ・実装演習(サポートベクターマシン)

はじめに

今回の記事は、前回に引き続き、私が受講しているE資格 JDLA認定プログラムの「ラビットチャレンジ」のレポート記事です。

今回のテーマは機械学習のSVM(サポートベクターマシン)です。
コード実装を取り入れながら行っていきます。

SVM(サポートベクターマシン)

回帰または分類に適用できる教師あり学習の機械学習手法です。
マージン最大化と呼ばれる考えに基づいて、主に2値の分類問題に使われます。
ここではその2値の分類問題について記載します。
分類の場合、線形モデルの正負で2値分類することになります。

線形判別関数(\(wTx+b=0wTx+b=0\))と最も近いデータ点(サポートベクトル)との距離(マージン)が最大となる線形判別関数を求めます。
平均二乗誤差(MSE)の残差は、グラフで表すと予測値からデータ点へ縦に真っ直ぐ差を取りますが、マージンは線形判別関数からデータ点へ90°の線を引いた差となります。

決定境界に近いデータの値が変わることにより、予測が間違うこともありますが、マージンを大きく取る(マージン最大化)によって、防ぐことができます。

マージン最大化

先述の通り、マージン最大化の目的は、マージン(直線とデータの距離)を一番大きくするベクトルを見つけることです。
この時の直線は、次のように表すことができます。
$$
W^TX_i + b
$$
マージンは、次のように表すことができます。
$$
M = \frac{W^TX_+ + b}{||W||} = \frac{-(W^TX_- + b)}{||W||}
$$
よって、マージン最大化は以下の最適化問題によって表すことができます。。
$$
max_{w,b} M,  \frac{t_i(W^TX_i + b)}{||W||},  i = 1,2,…,n
$$

線形分離可能と線形分離不可能

データを分類するための直線を引いた時、どのクラスに属するのかある程度分類できれば線形分離可能と言えます。
線形分離が可能な場合、ハードマージンを適用できます。
一方、どのクラスに属するのか不安定で、データにばらつきがある場合は線形分離不可能と言えます。
線形分離が不可能な場合、ソフトマージンが対応しています。
マージン最大化では、線形分離可能な場合と不可能な場合では計算式が異なります。

線形分離可能な場合

マージン最適化問題

$$
min_{w,b} \frac{1}{2}||W||^2\\
t_i(W^TX_i + b) ≤1\\
i = 1,2,…,n
$$

線形分離不可能な場合

マージン最適化問題

$$
min _{w, ξ} (\frac{1}{2}||W||^2 + C∑ξ_i)\\
\\
t_i(W^TX_i + b) ≧1ーξ_i ,ξ_i≧0
\\
i = 1,2,…,n
$$

双体表現

SVMのマージン最大化は、先述のハードマージン、ソフトマージンの最適化問題によって定式化できます。
これらの最適化問題を主問題と呼ぶ。
世の中には線形分離不可能な場合がたくさんあり、多くの場合SV分類の主問題を直接解くことはできません。

そこで、主問題を双体問題と呼ばれる別の形に表現して最適化を考えることで、非線形の分類が可能になります。
双体問題を理解するために必須なのが、ラグランジュ未定乗数法です。

ラグランジュ未定乗数法

制約付き最適化問題の代表的な解法です。
$$
L(x,a) = f(x) + ∑a_ig_i(x)
$$
SVMの主問題を双体問題にした後、ラグランジュの未定乗数法を用いて最適化問題を解きます。
$$
目的関数 min _{w, ξ} (\frac{1}{2}||W||^2 + C∑ξ_i)\\
\\
不等式制約 t_i(W^TX_i + b) ≧1ーξ_i   ,ξ_i≧0  i = 1,2,…,n
\\
$$

ラグランジュ関数(ラグランジュ乗数はα、βを用いる)

$$
L(w,b,ξ,a,β) = \frac{1}{2}||W||^2 + C∑ξ_i – ∑a_i(t_i(W^TX_i + b) – 1 + ξ_i) – ∑β_iξ_i
$$
これを用いると、元の最適化問題は、下記のような双体問題の式で書けます。
$$
max_a (L(a) = ∑a_i – \frac{1}{2}∑∑a_ia_jy_iy_jX^T_iX_j)\\
∑a_iy_i = 0 \\
0≦a_i≦C \\
i = 1,2,…,n
$$

カーネル法

カーネル法とは、低次元のデータを高次元に写像して分離する方法です。
これまで表現できた双体問題を今度はカーネル法を用いて解きます。
$$
k(x_i,x_j) = φ(X)^T_iφ(X)_j
$$
ここである特定の条件を満たす関数を用いるとφ(X)を直接計算せずに内積を計算することができる。
$$
φ(X)^T_iφ(X)_j = (X^Ty)^2
$$

SVM実装演習

必要なライブラリのインポート

まず、必要なライブラリをインポートしておきます。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

線形カーネル

線形分離可能なデータの生成

乱数に-2と+2を加算し、乱数を調整します。

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)

最急降下法(バッチ勾配降下法)

線形カーネルを使用します。

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)

予測

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.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')

マージンと決定境界を可視化すると、分離できていることがわかります。
丸がついている点がサポートベクトルです。

RBFカーネル(ガウシアンカーネル)

線形分離不可能なデータの生成

非線形分離のデータを生成します。

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)

RBFカーネルを用いた学習

元空間では線形分離できないので、特徴空間上で線形分離します。

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)

予測

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=['--', '-', '--'])

非線形で分離できていることがわかります。

ソフトマージンSVM

重なりのあるデータの生成

重なりがあるので、綺麗に分離できないデータです。

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)

マージン内部のデータや誤分類データにペナルティを与える最急降下法

ペナルティのパラメータ\(C\)は1にしています。
np.clip(a, 0, C)で制約条件をclipしています。

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)

予測

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=['--', '-', '--'])

これまでと異なり、マージン内にデータが入ることを許容して決定境界が引かれていることがわかります。