Support Vector Machine

Support Vector Machine#

マージン#

データがクラス$\{C_1, C_2\}$のどちらに含まれるかを判断する2クラス識別問題について考える。教師ラベルは$y \in \{+1, -1\}$であり、それぞれデータが$C_1, C_2$のどちらに含まれるかを示すとする。係数ベクトルを$w=(w_1, ..., w_d)^T$、バイアス項を$b$、特徴量ベクトルを$x=(x_1, ..., x_d)^T$とおくと、線形識別関数は

\[ f(x) = w^T x + b \]

と表すことができる。

識別境界（識別超平面）は$f(x)=0$となる位置に描かれるとし、クラス1を$f(x) >= 0$、クラス2を$f(x) < 0$で表現するように学習させるとする。例えば次の図のように、ある識別関数が存在したとする。

../_images/751fc07e6d5178c46dadb4a2106c4551cc2c8740823705810bf71ccd4173df33.png

訓練データ中に存在しなかったノイズがテストデータに含まれていた場合、ノイズの分だけ識別を誤りやすくなる。しかし、訓練データの点が識別超平面からある値$h > 0$よりも離れるように学習させれば、$h$より小さなノイズに対しては正しく識別できるようになる。

例えば、以下の図の(a)と(b)はいずれもサンプルをうまく分離できているものの、(a)よりも(b)のほうがデータ点と識別超平面の距離があり、ノイズに対してより頑健で望ましい分類器であると考えられる。

../_images/26fe24bf7b3892a42a1738b648e01dd59df9730972a1d2d7c8fd75aaf2e292c0.png

となれば、 「識別超平面が訓練データからもっとも離れるように（両クラスの中間になるように）学習させればよいのではないか」 という考えが湧く。

これがサポートベクターマシン（Support Vector Machine: SVM）の考え方である。

識別境界$f(x) = 0$と最も近い各クラスの訓練データの点を サポートベクトル （support vector）といい、サポートベクトルと識別境界との距離（識別境界と最も近いデータ点の距離）を マージン （margin）という。

ある識別関数に対してとれるマージンの大きさは、両クラスの学習データを識別関数の法線ベクトル上に射影した長さの最小値

\[ \rho(w) = \min_{x \in C_1} \frac{w^T x}{||w||} - \max_{x \in C_2} \frac{w^T x}{||w||} \]

の半分である。$\rho(w)$はクラス間マージンという。次の図中の2つの破線の間の距離が$\rho(w)$である

../_images/4f4601b1fe0fb610d43cbb1ba425a4aede425d111f155fd68f74225869020769.png

ハードマージンSVM#

学習データの集合を$\mathcal{D}_L = \{(y_i, x_i)\}(i=1,...,N)$とする。係数ベクトルはバイアス項$b$を外に出す形で、$w=(w_1, ..., w_d)^T$と表記する。特徴量ベクトルは$x=(x_1,...,x_d)^T$である。$y_i=\{-1, +1\}$は教師データで、学習データ$x_i\in \mathbb{R}^d$がどちらのクラスに属するかを示す。

線形識別関数のマージンを$\kappa$とすれば全ての学習データで

\[ |w^T x_i + b| \geq \kappa \]

が成り立つ。

係数ベクトルとバイアス項をマージンで正規化（$w^T x_i = -b$を定数倍）したものをあらためて$w, b$とおけば

\[\begin{split} \begin{cases} w^T x_i + b \geq +1 & \text{if} \hspace{0.5em} y_i = +1\\ w^T x_i + b \leq -1 & \text{if} \hspace{0.5em} y_i = -1 \end{cases} \end{split}\]

となり、まとめて表記すると

\[ y_i \times (w^T x_i + b) \geq 1 \]

クラス間マージンは $$ \rho(w, b) = \min_{x \in C_{y=+1}} \frac{w^T x}{||w||} - \max_{x\in C_{y=-1}} \frac{w^T x}{||w||} $$

第1項の分子は$w^T x_i + b \geq +1$の最小値が$w^T x_i + b = 1$であることから$\min w^T x_i = 1 - b$

第2項の分子は$w^T x_i + b \leq -1$の最大値が$w^T x_i + b = -1$であることから$\max w^T x_i = -1 - b$

であることを使えば

\[\begin{split} \begin{align} \rho(w, b) &= \min_{x \in C_{y=+1}} \frac{w^T x}{||w||} - \max_{x\in C_{y=-1}} \frac{w^T x}{||w||}\\ &= \frac{1 - b}{||w||} - \frac{-1 - b}{||w||}\\ &= \frac{1 + 1 - b + b}{||w||}\\ &= \frac{2}{||w||} \end{align} \end{split}\]

となる。

識別関数の最大マージンは最大クラス間マージンの半分であるため、$\frac{1}{||w||}$となる。

最適識別超平面#

最適な識別超平面は、「すべての訓練データを正しく識別できる」という制約条件

\[ y_i (w^T x_i + b) \geq 1 \hspace{1em} (i=1,...,N) \]

の下でマージン$\frac{1}{\|w\|}$を最大化した解として得られる。マージンの最大化は$\|w\|$の最小化と等しいため、

\[ w_0 = \min \|w\| \]

として求めることができる。これは次の不等式制約条件つき最適化問題を解くことで得られる。

主問題

\[\begin{split} \begin{align} \text{minimize} & \hspace{1em} L_p(w) = \frac{1}{2} w^T w\\ \text{subject to} & \hspace{1em} y_i( w^T x_i + b) \geq 1; \ \forall i \end{align} \end{split}\]

この問題はラグランジュの未定乗数法を用いて解かれ、次のラグランジュ関数として定式化される

\[ \tilde{L}_p(w, b, \alpha) = \frac{1}{2} w^T w - \sum^N_{i=1} \alpha_i (y_i (w^T x_i + b) - 1) \]

ここで$\alpha=(\alpha_1, ..., \alpha_N)^T$、$\alpha_i \geq 0$であり、$\alpha_i$はラグランジュ未定乗数と呼ばれる。

この最適化問題の解$w_*$と$b_*$は以下のKKT（Karush-Kuhn-Tucker）条件を満たす解として知られている。

KKT条件

(1) $\displaystyle \frac{\partial \tilde{L}_p(w, b, \alpha)}{\partial w}|_{w=w_*} = w_* - \sum^N_{i=1} \alpha_i y_i x_i = 0$

(2) $\displaystyle \frac{\partial \tilde{L}_p(w, b, \alpha)}{\partial b} = \sum^N_{i=1} \alpha_i y_i = 0$

(3) $y_i(w^T x_i + b) - 1 \geq 0$

(4) $\alpha_i \geq 0$

(5) $\alpha_i (y_i (w^T x_i + b) - 1) = 0$

ラグランジュ関数の$w$を$w_*$に置き換えてKKT条件(1)と(2)を代入して整理すると

\[\begin{split} \begin{align} L_d(\alpha) &= \frac{1}{2} {w_*}^T w_* - \sum^N_{i=1} \alpha_i y_i w_*^T x_i - b \sum^N_{i=1} \alpha_i y_i + \sum^N_{i=1} \alpha_i\\ &= \sum^N_{i=1} \alpha_i - \frac{1}{2} w_*^T w_* \hspace{2em} (\because \sum^N_{i=1} \alpha_i y_i = 0)\\ &= \sum^N_{i=1} \alpha_i - \frac{1}{2} \sum^N_{i=1} \sum^N_{j=1} \alpha_i \alpha_j y_i y_j x_i^T x_j \end{align} \end{split}\]

となり、ラグランジュ未定乗数のみの関数にすることができる

KKT条件(1)より最適解は$w_* = \sum^N_{i=1} \alpha_i y_i x_i$のようになることがわかっているので、最適な係数$\alpha_i$を求める問題に置き換えることができる。

双対問題

\[\begin{split} \begin{align} \text{maximize} & \hspace{1em} L_d(\alpha) = \alpha^T \boldsymbol{1} - \frac{1}{2} \alpha^T H \alpha \\ \text{subject to} & \hspace{1em} \alpha^T y = 0 \end{align} \end{split}\]

ここで

\[\begin{split} \begin{align} \boldsymbol{1} &= (1,...,1)^T\\ H &= (H_{ij} = y_i y_j x_i^T x_j)\\ y &= (y_1,...,y_N)^T\\ \end{align} \end{split}\]

である。

双対問題のラグランジュ関数は、ラグランジュ未定乗数を$\beta$とすれば次の関数になる。

\[ \tilde{L}_d(\alpha, \beta) = \alpha^T \boldsymbol{1} - \frac{1}{2} \alpha^T H \alpha - \beta \alpha^T y \]

KKT条件(5)より$\alpha_i (y_i (w^T x_i + b) - 1) = 0$がすべての$i$で成り立てば良いため、

\[\begin{split} \begin{cases} \alpha_i > 0 & \text{if} \hspace{0.5em} y_i(w^T x_i + b) - 1 = 0\\ \alpha_i = 0 & \text{if} \hspace{0.5em} y_i(w^T x_i + b) - 1 \neq 0 \end{cases} \end{split}\]

となる。$\alpha_i > 0$となる$x_i$をサポートベクトルという。

最適なバイアス$b_*$はサポートベクトルの一つ$x_s$を用いて

\[ y_s (w_*^T x_s + b_*) - 1 = 0 \]

を解いて求めるか、それらの平均を用いる。

実装例（cvxpy）#

主問題をそのままソルバーに通すパターン

\[\begin{split} \begin{align} \min_w & \hspace{1em} \frac{1}{2} w^T w = \frac{1}{2} \|w\|^2_2\\ \text{s.t.} & \hspace{1em} y_i( w^T x_i + b) \geq 1; \ \forall i \end{align} \end{split}\]

import cvxpy as cp
n = X.shape[0]
d = X.shape[1]

b = cp.Variable()
w = cp.Variable(d)
prob = cp.Problem(cp.Minimize( (1/2) * cp.norm(w, 2)**2 ),
                  [y[i] * (w.T @ X[i, :] + b) >= 1 for i in range(n)])
prob.solve()
print("The optimal value is", prob.value)
print("w is", w.value)
print("b is", b.value)

The optimal value is 0.9049773766503881
w is [0.90497738 0.99547511]
b is -8.920766380091388e-17

../_images/111df84fb84fbe8d4f70ecaff5c6a28d1125a0e9dcbbc4e6dd267af1b19c46c4.png

双対問題をソルバーに通すパターン#

双対問題

\[\begin{split} \begin{align} \text{maximize} & \hspace{1em} L_d(\alpha) = \alpha^T \boldsymbol{1} - \frac{1}{2} \alpha^T H \alpha \\ \text{subject to} & \hspace{1em} \alpha^T y = 0 \end{align} \end{split}\]

をcvxpyの二次計画問題のソルバーを使って解いてみる

# データ数が多くなるとソルバーが上手く動かないので一旦暫定的対処としてデータ数を絞る
y = y[3:6]
X = X[3:6]

# Hを作成
n = X.shape[0]
H = np.zeros(shape=(n, n))
for i in range(n):
    for j in range(n):
        H[i, j] = y[i] * y[i] * X[i] @ X[j]

ones = np.ones(shape=(n, ))

import cvxpy as cp

alpha = cp.Variable(n)  # 長さnのベクトル
prob = cp.Problem(cp.Maximize( alpha.T @ ones - (1/2) * cp.quad_form(alpha, H) ),
                  [alpha.T @ y == 0])
prob.solve()
print("The optimal value is", prob.value)
print("alpha is", alpha.value)

a = alpha.value
w = sum([a[i] * y[i] * X[i] for i in range(n)])
print(f"w={w}")

The optimal value is 33.99999999999994
alpha is [34. 10. 24.]
w=[39.6 49. ]

../_images/e7b6f57c1692c8615c0cc2be8f8519f32e41cf0f1309bb0c4601fb58c5c98a2b.png

実装例（scikit-learn）#

from sklearn import svm
clf = svm.SVC(random_state=0, kernel='linear', shrinking=False)
clf.fit(X, y)

print(f"b={clf.intercept_}, w={clf.coef_}")

b=[-0.15974441], w=[[0.76677316 0.83067093]]

../_images/92b9dbca1faf98f0deb48199e7f905176dc9bbe1c74718b1c21d06efe6bae0a6.png

ソフトマージンSVM#

C-SVM#

スラック変数と呼ばれる変数$\xi_i$を追加する。

\[\begin{split} \begin{cases} \xi_i = 0 & (マージン内で正しく識別できる場合)\\ 0 < \xi_i \leq 1 & (マージン境界を超えるが正しく識別できる場合)\\ \xi_i > 1 & (識別境界を超えて誤識別される場合) \end{cases} \end{split}\]

以下のように書くこともできる

\[ \xi_i = \max[0, 1-y_i(w^T x_i + b)] = f_{+}(1 - y_i(w^T x_i + b)) \]

ここで$f_{+}(x)$はヒンジ（hinge）関数と呼ばれるもので

\[\begin{split} f_{+}(x) := \begin{cases} x & (x > 0の場合)\\ 0 & (それ以外) \end{cases} \end{split}\]

である

ソフトマージン識別器の主問題は以下のように定式化される。

主問題

\[\begin{split} \begin{align} \text{minimize} & \hspace{1em} L_p(w, \xi) = \frac{1}{2} w^T w + C\sum^N_{i=1} \xi_i \\ \text{subject to} & \hspace{1em} y_i( w^T x_i + b) - 1 + \xi_i \geq 0 \\ & \hspace{1em} \xi_i \geq 0 \end{align} \end{split}\]

すべての訓練データのスラック変数の和$\sum \xi_i (\xi_i \geq 0)$は誤識別数の上限を与える。パラメータ$C$は誤識別数に対するペナルティの強さであり、$C$が大きいほど$w$のノルム最小化よりも誤識別数を小さくする方を優先することになる。

このSVMは$C$-SVMと呼ばれる。

ν-SVM#

上限サポートベクトル（マージン誤り$\xi_i > 0$のベクトルの数）の割合の上限を規定するハイパーパラメータ$\nu$が指定できるようになった

カーネルトリック#

カーネルモデル#

線形モデルをカーネルモデルに拡張することを考える。

線形モデル

\[ f^{linear}(x_i) = w^T x_i + b \]

$x_i = (x_{i1}, \dots, x_{iD})$
$w=(x_1,\dots,x_D)^T$

カーネルモデル

\[ f^{kernel}(x_i) = \sum^N_{j=1} w_j \phi(x_i)^T \phi(x_j) + b \]

$x_i = (x_{i1}, \dots, x_{iD})$
$w=(x_1,\dots,x_N)^T$

ここで$\phi(\cdot)$は任意の関数で、$\phi(\cdot)$によって入力ベクトル$x$を高次元空間に写像し（2次元では線形分離不可能なものを3次元に写して線形分離不可にするイメージ）、高次元空間上の類似度を内積$\phi(x_i)^T \phi(x)$で表す。

カーネルモデルは$N$の和が入っているように、訓練データ数$N$が増えるとモデルの表現力は高まるが計算量が増える。

カーネルトリック#

$\phi(x)$上での内積計算を緩和するために 正定値関数 （positive definite function）を用いる。

正定値関数

関数 $k(x_i, x_j)$ は次の条件を満たすとき正定値関数と呼ばれる。

(1) 対称性：$k(x_i, x_j)=k(x_j, x_i)$

(2) 正定値性：デー夕点 $\mathrm{x}_1, \mathrm{x}_2, \ldots, \mathrm{x}_N$ に対する以下の グラム行列（Gram matrix） $K$ が半正定値である。

\[\begin{split} K = \begin{pmatrix} k(x_1, x_1) & k(x_1, x_2) & \cdots & k(x_1, x_N) \\ k(x_2, x_1) & k(x_2, x_2) & \cdots & k(x_2, x_N) \\ \vdots & \vdots & \ddots & \vdots \\ k(x_N, x_1) & k(x_N, x_2) & \cdots & k(x_N, x_N) \end{pmatrix} \end{split}\]

すなわち、任意の $N$ 次元のベクトル $z$ に対し $z^{\top} K z \geq 0$ が成り立つ。

正定値関数 $k(x_i, x_j)$ は 再生核ヒルベルト空間 $\mathcal{H}_k$への写像 $\phi(x_i), \phi(x_j) \in \mathcal{H}_k$の内積に対応する。

\[ k(x_i, x_j) = \phi(x_i)^T \phi(x_j) \]

これを用いて高次元空間上での内積をより単純な2変数関数の計算$k(\cdot, \cdot)$に置き換える事ができる。この正定値関数$k(\cdot, \cdot)$のことを カーネル関数 （kernel function）、グラム行列$K$を カーネル行列 （kernel matrix）と呼ぶ。

カーネル関数とカーネル行列を用いると、カーネルモデルは以下のように表現できる

カーネル行列を用いたカーネルモデル

\[ f^{kernel}(x_i) = \sum^N_{j=1} w_j k(x_i, x_j) + b = K_{i:} w + b \]

ここで$K_{i:}$はカーネル行列の$i$行目の行ベクトルを表す。

カーネル関数の例#

線形カーネル（linear kernel）

\[ k(x_i, x_j) = x_i x_j^T \]

入力をそのまま出力する写像関数$\phi(x)=x^T$に対応する

ガウスカーネル（Gaussian kernel）

\[ k(x_i, x_j) = \exp \left( -\frac{\| x_i-x_j \|^2}{2\sigma^2} \right) \]

カーネルモデルを用いたソフトマージン最大化問題（主問題）

\[\begin{split} \begin{align} \min_{w,b,\xi} & \hspace{1em} L_p(w, \xi) = \frac{1}{2} w^T w + C\sum^N_{i=1} \xi_i \\ \text{s. t.} & \hspace{1em} y_i ( K_{i:} w + b) - 1 + \xi_i \geq 0 \\ & \hspace{1em} \xi_i \geq 0 ~ \forall i \end{align} \end{split}\]

カーネルモデルを用いたソフトマージン最大化のラグランジュ双対問題

\[\begin{split} \begin{aligned} \max _{\boldsymbol{\lambda}} \mathcal{L}\left(\boldsymbol{\lambda}, \boldsymbol{\gamma}^*, \mathbf{w}^*, b^*, \boldsymbol{\xi}^*\right) = & -\frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N \lambda_i y_i k (x_i, x_j) t_j \lambda_j+\sum_{i=1}^N \lambda_i \\ \text { s.t. } \quad & -\lambda_i \leq 0 \\ & \lambda_i \leq C \\ & \sum_{i=1}^N \lambda_i y_i=0 \forall_i \end{aligned} \end{split}\]

# TODO: コードかく

参考#

八谷大岳. (2020). ゼロからつくるPython機械学習プログラミング入門.

# TODO

from sklearn import datasets
X, y = datasets.load_iris(return_X_y=True, as_frame=False)
X = X[:, 0:2]

fig, ax = plt.subplots()
for c, label in enumerate(['setosa', 'versicolor', 'virginica']):
    is_c = y == c
    ax.scatter(X[is_c, 0], X[is_c, 1], label=label)
ax.legend()

<matplotlib.legend.Legend at 0x7f4cca912d10>

../_images/a1be46ec7418da0936783b58f7aeb9ecde8860397b1d618e9b132dbd17fd70fa.png

# setosaとvirginicaを1つにまとめて2クラスにする
y[y == 2] = 0


fig, ax = plt.subplots()
for c, label in enumerate(['setosa', 'versicolor', 'virginica']):
    is_c = y == c
    ax.scatter(X[is_c, 0], X[is_c, 1], label=label)
ax.legend()

<matplotlib.legend.Legend at 0x7f4cd4ba76d0>

../_images/d1f3b5b14de9c59c0c3811a5c2b95f41b02f4d70832c79e5a035410a0ed3c7da.png

GaussianKernel (RBF kernel)#

一般的には $$ k(\mathbf{x},\mathbf{x}^{\prime}) = \exp \left( \frac{-\| \mathbf{x} - \mathbf{x}^{\prime} \|^{2}}{2\sigma^{2}} \right) \tag{2} $$

David Duvenaud (2014). “The Kernel Cookbook: Advice on Covariance functions”.

x1 = np.array([[1], [2]])
x2 = np.array([[3], [4]])
x1

array([[1],
       [2]])

def rbf_kernel(x1, x2, sigma=1.0) -> float:
    return - np.linalg.norm(x1 - x2)**2 / (2 * sigma**2)
rbf_kernel(x1, x2)

-4.000000000000001

class GaussianKernel:

    def fit(self, X):
        self.alpha = 0.001
        self.mu = np.mean(X, axis=0)
        self.sigma = np.cov(X, rowvar=False, bias=False)

    def transform(self, X) -> np.array:
        return np.array([self._transform(x) for x in X])

    def _transform(self, x) -> float:
        return np.exp( - self.alpha * (x - self.mu) @ np.linalg.inv(self.sigma) @ (x - self.mu) )
    
    def __repr__(self) -> str:
        return f"<{self.__class__.__name__} alpha={self.alpha}, mu={self.mu}, sigma={self.sigma}>"

Memo 4#

RBFカーネルが無限次元になるのは指数関数の冪級数による定義が無限和であるため

\[ \exp(x) = \sum^\infty_{n=0} \frac{1}{n!} x^n \]

ref: https://ja.wikipedia.org/wiki/指数関数#厳密な定義a

TODO: もっと図とか入れたい

9. サポートベクトルマシン — 機械学習帳

のように

Support vector machine classifier with (\ell_1)-regularization — CVXPY 1.3 documentation

はじパタ全力解説: 第8章サポートベクトルマシン - Qiita

Support Vector Machine

Contents

Support Vector Machine#

マージン#

ハードマージンSVM#

最適識別超平面#

実装例（cvxpy）#

双対問題をソルバーに通すパターン#

実装例（scikit-learn）#

ソフトマージンSVM#

C-SVM#

ν-SVM#

カーネルトリック#

カーネルモデル#

カーネルトリック#

カーネル関数の例#

参考#

GaussianKernel (RBF kernel)#

Memo 4#