In [1]:

import numpy as np 
import matplotlib.pyplot as plt

In [2]:

seed = 0

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

マージン最大化に基づく線形識別器
カーネル関数と組み合わせることによって非線形な識別も可能
サンプル間の類似度にしか依存しない（→非ベクトルデータの分類にも適用可能）

まとめ

マージン最大化という方針で誤差関数を定義する。
汎化性能を高めるために、スラック変数を導入してソフトなマージン最大化に緩和する
ラグランジュ未定乗数法を用いてラグランジュ関数を得る。
各パラメータについて微分する。
ラグランジュ関数に代入して双対問題を得る。
SMO法などを用いて二次計画問題である双対問題を解く。
無効制約のデータ点は識別関数に影響しなくなる。（SVMが疎になる。）

マージン最大化¶

「個々のデータに関して確率計算をし、何らかの識別面を求める」のではなく、「汎化性能の最も高い識別面を直接考える」
そこで、「マージン最大化」という方針で識別面を最適化する。（正則化理論・VC理論により正当化されている。）
最も近いデータ $\phi(\mathbf{x}_n)$ までの距離（マージン）が最大となるような識別面を選ぶ。（ここでは、距離＝ユークリッド距離とする。）
識別面からちょうどマージンの分だけ離れた面上に乗る特徴ベクトル（一般に複数）を、(広義の) サポートベクトル (support vector) と呼ぶ。

学習法の定式化（2クラス分類）¶

※ 二クラスが線形分離可能であると仮定する。

線形識別関数が $y(\mathbf{x}) = \mathbf{w}^T\phi(\mathbf{x})+b$ で表され、$y(\mathbf{x})$ の正負でクラスを識別するとする。
- $y(\mathbf{x})>0$ の時 $t=1$
- $y(\mathbf{x})<0$ の時 $t=-1$
したがって、$t_ny(\mathbf{x}_n)>0$ が全てのデータ $(\mathbf{x}_n,t_n)$ に対して成り立つ（∵特徴空間で線形分離可能）
学習の目的は「パラメータ $\mathbf{w}$ と $b$ を決定すること」。
データと識別面の最小距離は $1 / \|\mathbf{w}\|$ で表される。（理由は以下）

※ 記述を簡単にするために線形識別問題の場合で考えるが、非線形写像 $\phi$ を通しても同様の議論が成り立つ。

点 $\mathbf{x}_n$ から識別面に下ろした垂線の足を $\mathbf{x}_{n\perp}$
点 $\mathbf{x}_n$ から識別面までの符号付き距離を $d_n$

と定義すると、

$\mathbf{w}$ が識別面の法ベクトルであるから、 $\mathbf{x}_{n\perp} = \mathbf{x}_n - d_n\frac{\mathbf{w}}{\|\mathbf{w}\|}$ と表せる。
$\mathbf{x}_{n\perp}$ は識別面上にあるので、$\mathbf{w}^T\mathbf{x}_{n\perp} + b = 0$ が成り立つ。

したがって、 $$\mathbf{w}^T\left(\mathbf{x}_n - d_n\frac{\mathbf{w}}{\|\mathbf{w}\|}\right) + b = 0\qquad \therefore d_n = \frac{\mathbf{w}^T\mathbf{x}_n + b}{\|\mathbf{w}\|}$$ として距離 $d_n$ が定義できる。ここで、$(\mathbf{w},b)\rightarrow(k\mathbf{w},kb)$ としても $d_n$ の値は不変なので、

$$t_n\left(\mathbf{w}^T\mathbf{x}_n + b \right)= 1\qquad (7.4)$$

とスケーリングしても一般性を失わない。なお、このスケーリング下では、全てのデータ $\mathbf{x}_n$ に対して

$$t_n\left(\mathbf{w}^T\mathbf{x}_n + b \right) \geq 1\qquad (7.5)$$

が成り立つ。</span>

ここまでをまとめると、マージン最大は、『$t_n\left(\mathbf{w}^T\mathbf{x}_n+b\right)\geq1\quad(i=1,\ldots,n)$ という条件下で、$\frac{1}{2}||\mathbf{w}||^2$ を最小化する問題』であると言える。

この制約付き最適化問題を解くために、$(7.5)$ の各制約式ごとにラグランジュ乗数 $a_n\geq0$ を導入すると、次のラグランジュ関数が得られる。

$$L(\mathbf{w},b,\mathbf{a}) = \frac{1}{2}\|\mathbf{w}\|^2 - \sum_{n=1}^{N}a_n\left\{t_n\left(\mathbf{w}^T\phi(\mathbf{x}_n) + b\right) - 1\right\}\qquad (7.7)$$

よって、$\mathbf{w},b$ について微分して、

$$ \begin{aligned} \frac{\partial L}{\partial \mathbf{w}} &= \mathbf{w}-\sum_{n=1}^N a_nt_n\phi(\mathbf{x}_n) =0 \quad \therefore \mathbf{w}^{\ast}=\sum_{n=1}^N a_nt_n\phi(\mathbf{x}_n) & (7.8)\\ \frac{\partial L}{\partial b} &= -\sum_{n=1}^N a_nt_n = 0 \quad\therefore\sum_{n=1}^Na_n^{\ast}t_n = 0 & (7.9)\\ \end{aligned} $$

が得られる。これらの結果を $(7.7)$ に代入して $\mathbf{w},b$ を消去すると、元の問題は以下の双対問題に変換される。

$$ \tilde{L}(\mathbf{a}) =\sum_n^Na_n - \frac{1}{2}\sum_{n=1}^{N}\sum_{m=1}^Na_na_mt_nt_mk(\mathbf{x}_n,\mathbf{x}_m)\qquad (7.10)$$

なお、$\mathbf{a}$ は以下の制約を満たす。

$$ \begin{aligned} a_n & \geq 0,\quad n=1,\ldots,N & (7.11)\\ \sum_{n=1}^Na_nt_n &= 0 & (7.12) \end{aligned} $$

また、カーネル関数 $k(\mathbf{x},\mathbf{x}^{\prime})$ は、$k(\mathbf{x},\mathbf{x}^{\prime}) = \phi(\mathbf{x})^T\phi(\mathbf{x})$ と定義される。

この双対問題は二次計画問題であり、逐次最小問題最適化法(SMO;sequential minimal optimization)などで解くことができる。

なお、今回の場合、制約条件が全て線形であることから実行可能域が凸領域であることがわかるので、任意の局所解が大局解になる。

一旦最適な $\mathbf{a}$ が求まると、識別関数 $y(\mathbf{x})$ は、

$$y(\mathbf{x}) = \mathbf{w}^{\ast T}\phi(\mathbf{x})+b = \sum_{n=1}^N a_nt_nk(\mathbf{x},\mathbf{x}_n) + b \qquad (7.13)$$

として表すことができる。ここで、KKT条件より $a_n=0\Leftrightarrow t_ny(\mathbf{x}_n)-1>0$ のデータは識別関数に影響しないので、計算から除外することができる。（SVMが疎になる理由）

なお、バイアスパタメータ $b$ は、任意のサポートベクトル $\mathbf{x}_n$ が $t_ny(\mathbf{x}_n) = 1$ を満たすことから、

$$ \begin{aligned} &t_n\left(\sum_{m\in\mathcal{S}}a_mt_mk(\mathbf{x}_n,\mathbf{x}_m) + b\right) = 1 & (7.17)\\ &\therefore b^{\ast}= \frac{1}{N_{\mathcal{S}}}\sum_{n\in\mathcal{S}}\left(t_n-\sum_{n\in\mathcal{S}}a_mt_mk(\mathbf{x}_n,\mathbf{x}_m)\right)&(7.18) \end{aligned} $$

として求めることができる。（数値計算の誤差の影響を小さくするため、$(7.17)$ の両辺に $t_n$ をかけ、全てのサポートベクトルで平均をとった。

In [3]:

from kerasy.ML.svm import hardSVC

In [4]:

n_cls1=20; n_cls2=20
N = n_cls1+n_cls2

In [5]:

train_x = np.concatenate([
    np.random.RandomState(seed).normal(loc=[0.5,  0], scale=[0.2,0.4], size=(n_cls1,2)),
    np.random.RandomState(seed).normal(loc=[-0.5, 0], scale=[0.2,0.4], size=(n_cls2,2))
])
train_y = np.concatenate([
    np.ones(shape=n_cls1),
    -np.ones(shape=n_cls2)
])

In [6]:

X,Y = np.meshgrid(np.linspace(-1, 1, 100), np.linspace(-1, 1, 100))

In [7]:

kernels = [
    "linear","polynomial","gaussian",
    "sigmoid","laplacian","inverse_multiquadric"
]

In [8]:

n_fig = len(kernels)
n_col = 3
n_row = n_fig//n_col if n_fig%n_col==0 else n_fig//n_col+1

fig = plt.figure(figsize=(4*n_col,4*n_row))
for i,kernel in enumerate(kernels):
    ax = fig.add_subplot(n_row,n_col,i+1)
    
    model = hardSVC(kernel=kernel)
    model.fit(train_x, train_y, max_iter=1500, sparse_memorize=False)
    
    Z = np.vectorize(lambda x,y: model.predict(np.array([x,y]).reshape(1,-1)))(X, Y) == 1
    ax.pcolor(X, Y, Z, alpha=0.3)
    for i in range(N):
        alpha = 1 if i in model.isSV() else 0.3
        color = "red" if train_y[i]==1 else "blue"
        ax.plot(train_x[i,0], train_x[i,1], marker="o", color=color, markersize=10, alpha=alpha)
    ax.set_title(f"{kernel} kernel")
    ax.set_xlim(-1,1), ax.set_ylim(-1,1)
plt.tight_layout()
plt.show()

ソフトマージン¶

ここまで、「二クラスが線形分離可能である」という仮定を置いていたが、実際の問題ではクラスの分布が重なる場合があり、そのような場合訓練データを完全に分離する解が汎化能力に優れるとは限らない（というより、過学習する。）

そこで、各学習データ $(\mathbf{x}_n,t_n)$ に対して、スラック変数 (slack variable) $\xi_n\geq 0$ を導入して、

$$ t_ny(\mathbf{x}_n) \geq 1 \Longrightarrow t_ny(\mathbf{x}_n) \geq 1 - \xi_n \qquad (7.20)$$

と条件を緩める。（ハードマージンの制約のソフトマージンへの緩和）つまり、$\xi_n > 1$ となるデータは誤識別している事になる。

ここでの目標は『マージンよりも内側（$\xi_n>0$）のデータ $n$ に対してソフトにペナルティを与えつつマージンを最大化する』ことである。したがって、次の目的関数を最小化する。

$$C\sum_{n=1}^N\xi_n + \frac{1}{2}\|\mathbf{w}\| \qquad (7.21)$$

なお、$C>0$ はスラック変数を用いて表されるペナルティとマージンの大きさの間のトレードオフを制御するパラメータである。

したがって、この最小化問題のラグランジュ関数は、新たにラグランジュ乗数 $\{\mu_n\}$ を導入して以下のようになる。

$$L(\mathbf{w},b,\boldsymbol{\xi},\mathbf{a},\boldsymbol{\mu}) = \frac{1}{2}\|\mathbf{w}\|^2 + C\sum_{n=1}^N\xi_n - \sum_{n=1}^{N}a_n\left\{t_n\left(\mathbf{w}^T\phi(\mathbf{x}_n) + b\right) - 1 + \xi_n\right\} - \sum_{n=1}^N\mu_n\xi_n\qquad (7.22)$$

先程同様 $\mathbf{w},b,\xi_n$ について微分して、

$$ \begin{aligned} \frac{\partial L}{\partial \mathbf{w}} &= \mathbf{w}-\sum_{n=1}^N a_nt_n\phi(\mathbf{x}_n) =0 \quad \therefore \mathbf{w}^{\ast}=\sum_{n=1}^N a_nt_n\phi(\mathbf{x}_n) & (7.29)\\ \frac{\partial L}{\partial b} &= -\sum_{n=1}^N a_nt_n = 0 \quad\therefore\sum_{n=1}^Na_n^{\ast}t_n = 0 & (7.30)\\ \frac{\partial L}{\partial \xi_n} &= C-a_n-\mu_n \quad\therefore a_n = C-\mu_n & (7.31) \end{aligned} $$

が得られる。これらの結果を $(7.22)$ に代入すると、以下の双対形のラグランジュ関数が得られる。

$$ \tilde{L}(\mathbf{a}) =\sum_n^Na_n - \frac{1}{2}\sum_{n=1}^{N}\sum_{m=1}^Na_na_mt_nt_mk(\mathbf{x}_n,\mathbf{x}_m)\qquad (7.32)$$

これは、先程の $(7.10)$ と同様の形である。なお、制約式は

$$ \begin{aligned} &a_n \geq 0,\mu_n\geq 0 \quad n=1,\ldots,N\\ &\therefore 0\leq a_n \leq C & (7.33)\\ &\sum_{n=1}^Na_nt_n = 0 & (7.34) \end{aligned} $$

である。

なお、バイアスパタメータ $b$ は、サポートベクトルは $a_n>0$ を満たすが、KKT条件より

$a_n$	$\mu_n$	$\xi_n$	説明
$a_n=0$			識別関数に影響を及ぼさない。
$0 < a_n < C$	$\mu_n>0$	$\xi_n=0$	ちょうどマージン境界上ある。
$a_n=C$	$\mu_n=0$	$\xi_n>0$	マージン内に侵入している。

となることから、$0 < a_n < C$ となるサポートベクトルについては $t_ny(\mathbf{x}_n) = 1$ が成立するので、

$$ \begin{aligned} &t_n\left(\sum_{m\in\mathcal{S}}a_mt_mk(\mathbf{x}_n,\mathbf{x}_m) + b\right) = 1 & (7.36)\\ &\therefore b^{\ast}= \frac{1}{N_{\mathcal{M}}}\sum_{n\in\mathcal{M}}\left(t_n-\sum_{n\in\mathcal{S}}a_mt_mk(\mathbf{x}_n,\mathbf{x}_m)\right)&(7.37) \end{aligned} $$

として求めることができる。

In [9]:

from kerasy.ML.svm import SVC

In [10]:

train_x = np.concatenate([
    np.random.RandomState(seed).normal(loc=[0.2,  0.2], scale=[0.2,0.4], size=(n_cls1,2)),
    np.random.RandomState(seed).normal(loc=[-0.2, -0.2], scale=[0.3,0.2], size=(n_cls2,2))
])
train_y = np.concatenate([
    np.ones(shape=n_cls1),
    -np.ones(shape=n_cls2)
])

In [11]:

params = [
    (1e-2,1e-2),(1e-2,1e0),(1e-2,1e2),
    (1e1, 1e-2),(1e1, 1e0),(1e1, 1e2),
    (1e4, 1e-2),(1e4, 1e0),(1e4, 1e2),
]

In [12]:

n_fig = len(params)
n_col = 3
n_row = n_fig//n_col if n_fig%n_col==0 else n_fig//n_col+1

fig = plt.figure(figsize=(4*n_col,4*n_row))
for i,(C,sigma) in enumerate(params):
    ax = fig.add_subplot(n_row,n_col,i+1)
    
    model = SVC(kernel="gaussian",C=C,sigma=sigma)
    model.fit(train_x, train_y, max_iter=1500, sparse_memorize=False)
    
    Z = np.vectorize(lambda x,y: model.predict(np.array([x,y]).reshape(1,-1)))(X, Y) == 1
    ax.pcolor(X, Y, Z, alpha=0.3)
    for i in range(N):
        alpha = 1 if i in model.isSV() else 0.3
        color = "red" if train_y[i]==1 else "blue"
        ax.plot(train_x[i,0], train_x[i,1], marker="o", color=color, markersize=10, alpha=alpha)
    ax.set_title(f"gaussian kernel (C={C}, $\sigma={sigma}$)")
    ax.set_xlim(-1,1), ax.set_ylim(-1,1)
plt.tight_layout()
plt.show()

多クラスSVM¶

SVMは本来2クラス識別の為の分類器であるが、これを複数組み合わせるなどの方法によって $K>2$クラスの識別に利用する事が出来る。最もよく用いられる方法は、one-versus-the-rest方式と呼ばれる方法であり、「データ $\mathbf{x}$ がクラス $k$ かそれ以外であるかを識別するSVM $y_k(\mathbf{x})$ を $K$ クラスそれぞれについて構築し、$ y(\mathbf{x}) = \mathop{\rm arg~max}\limits_{k}y_k(\mathbf{x})$ によって $\mathbf{x}$ のクラスを決定する」

SVM回帰¶

単純な線形回帰問題においては、次の式で定義される正則化された誤差関数を最小化する。

$$\frac{1}{2}\sum_{n=1}^N\{y_n-t_n\}^2 + \frac{\lambda}{2}\|\mathbf{w}\|^2\qquad (7.50)$$

疎な解を得るために、二条誤差関数を $\epsilon$ 許容誤差関数($\epsilon$-insentive error function)で置き換える。つまり、予測値 $y(\mathbf{x})$ と観測値 $t$ の差が $\epsilon(>0)$ 未満のときは、$\epsilon$ 許容誤差関数の値は $0$ となる。$\epsilon$ 許容誤差関数の簡単な例は、「許容範囲外の値に関しては線形のコストを生じさせる」ものである。

$$ \begin{aligned} E_{\epsilon}\left(y(\mathbf{x})-t\right) = \begin{cases} 0 & (\text{ when }|y(\mathbf{x} - t)| < \epsilon)\\ |y(\mathbf{x} - t)| - \epsilon & (\text{ otherwise }) \end{cases}\qquad (7.51) \end{aligned} $$

In [13]:

squareError  = lambda x,y_true: 1/2*(x-y_true)**2 # Sum of Squares Error function
epsilonError = lambda x,y_true,epsilon=1: 0 if abs(x-y_true)<epsilon else abs(x-y_true)-epsilon

In [14]:

X = np.linspace(-2,2,1000)
plt.plot(X,[squareError(x,y_true=0)  for x in X], color="black", label="sum of squares error function", linestyle="--")
plt.plot(X,[epsilonError(x,y_true=0) for x in X], color="black", label="$\epsilon$-insentive error function")
plt.scatter(0,0,marker="x",color="red",label="y_true")
plt.title("Error function comparison", fontsize=16), plt.xlim(-2,2), plt.ylim(-.1,2)
plt.xlabel("$x$", fontsize=16), plt.ylabel("$E(x)$", fontsize=16)
plt.legend()
plt.show()

この $\epsilon$ 許容誤差関数を用いると、結局次の正則化された誤差関数を最小化することになる。

$$C\sum_{n=1}^NE_{\epsilon}\left(y(\mathbf{x}_n)-t_n\right) + \frac{1}{2}\|\mathbf{w}\|^2\qquad (7.52)$$

分類問題と同様、この最適化問題をスラック変数 $\xi_n,\hat{\xi}_n$ を導入して表現する。

$\epsilon$ チューブの中に位置するという条件は $y_n-\epsilon\leq t_n \leq y_n+\epsilon$ と表せるので、チューブの外側にデータ点が存在することを許すような制約条件は以下のようにかける。

$$ \begin{aligned} t_n&\leq y(\mathbf{x}_n) + \epsilon + \xi_n & (7.53)\\ t_n&\geq y(\mathbf{x}_n) - \epsilon - \xi_n & (7.54)\\ \end{aligned} $$

よって、スラック変数を用いると、SVM回帰の誤差関数は次のようにかける。

$$C\sum_{n=1}^N\left(\xi_n + \hat{\xi}_n\right) + \frac{1}{2}\|\mathbf{w}\|^2\qquad (7.55)$$

したがって、この最小化問題のラグランジュ関数は、ラグランジュ乗数 $a_n,\hat{a}_n,\mu_n,\hat{\mu}_n$ を導入して以下のようになる。（$y(\mathbf{x}_n)= \mathbf{w}^T\phi(\mathbf{x}_n)+b$）

$$ \begin{aligned} L &=C\sum_{n=1}^N\left(\xi_n+\hat{\xi}_n\right)+\frac{1}{2}\|\mathbf{w}\|^2-\sum_{n=1}^N\left(\mu_n\xi_n+\hat{\mu}_n\hat{\xi}_n\right)\\ &-\sum_{n=1}^Na_n\left(\epsilon + \xi_n + y_n - t_n \right)-\sum_{n=1}^N\hat{a}_n\left(\epsilon + \hat{\xi}_n - y_n + t_n\right) \qquad (7.56)\\ \end{aligned} $$

これまで同様 $\mathbf{w},b,\xi_n,\hat{\xi}_n$ について微分して、

$$ \begin{aligned} \frac{\partial L}{\partial \mathbf{w}} &= \mathbf{w}-\sum_{n=1}^Na_n\phi(\mathbf{x}_n)-\sum_{n=1}^N\hat{a}_n\phi(\mathbf{x}_n)=0 \quad \therefore \mathbf{w}^{\ast}=\sum_{n=1}^N \left(a_n-\hat{a}_n\right)\phi(\mathbf{x}_n) & (7.57)\\ \frac{\partial L}{\partial b} &= -\sum_{n=1}^N a_n + \sum_{n=1}^N\hat{a}_n = 0 \quad\therefore\sum_{n=1}^N\left(a_n-\hat{a}_n\right) = 0 & (7.58)\\ \frac{\partial L}{\partial \xi_n} &= C-a_n-\mu_n \quad\therefore a_n + \mu_n = C & (7.59)\\ \frac{\partial L}{\partial \hat{\xi}_n} &= C-\hat{a}_n-\hat{\mu}_n \quad\therefore \hat{a}_n + \hat{\mu}_n = C & (7.60)\\ \end{aligned} $$

が得られる。これらの結果を $(7.56)$ に代入すると、以下の双対形のラグランジュ関数が得られる。

$$ \begin{aligned} \tilde{L}(\mathbf{a},\hat{\mathbf{a}}) &= -\frac{1}{2}\sum_{n=1}^N\sum_{m=1}^N\left(a_n-\hat{a}_n\right)\left(a_m-\hat{a}_m\right)k(\mathbf{x}_n,\mathbf{x}_m)\\ &- \epsilon\sum_{n=1}^N\left(a_n + \hat{a}_n\right) + \sum_{n=1}^N\left(a_n-\hat{a}_n\right)t_n \end{aligned} \qquad (7.61) $$

なお、制約式は

$$ \begin{aligned} 0 \leq a_n &\leq C & (7.62)\\ 0 \leq \hat{a}_n &\leq C & (7.63) \end{aligned} $$

である。したがって、新しい入力に対する予測値は

$$y(\mathbf{x}) = \sum_{n=1}^N\left(a_n-\hat{a}_n\right)k(\mathbf{x},\mathbf{x}_n) + b\qquad (7.64)$$

とカーネル関数のみを用いて計算できることがわかる。

なお、バイアスパタメータ $b$ は次のように計算できる。

$C>a_n>0$ が成り立つデータ点についてはKKT条件より $\xi_n = 0 \quad \therefore \epsilon+y_n-t_n=0$ が成り立つ。よって、 $$ \begin{aligned} b &= t_n - \epsilon - \mathbf{w}^T\phi(\mathbf{x}_n)\\ &= t_n - \epsilon - \sum_{m=1}^N\left(a_m-\hat{a}_m\right)k(\mathbf{x}_n,\mathbf{x}_m)\qquad (7.69) \end{aligned} $$
$C>\hat{a}_n>0$ が成り立つデータ点についてはKKT条件より $\hat{\xi}_n = 0 \quad \therefore \epsilon+y_n-t_n=0$ が成り立つ。よって、 $$ \begin{aligned} b &= t_n - \epsilon - \mathbf{w}^T\phi(\mathbf{x}_n)\\ &= t_n - \epsilon - \sum_{m=1}^N\left(a_m-\hat{a}_m\right)k(\mathbf{x}_n,\mathbf{x}_m)\qquad (7.69) \end{aligned} $$

なお、このようにして得られた $b$ の値を全て平均した結果を用いると数値誤差の影響を小さくできる。

In [ ]:

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