Support Vector Machine

Support Vector Machine#

マージン#

データがクラス$\{C_1, C_2\}$のどちらに含まれるかを判断する2クラス識別問題について考える。

教師ラベルは$t \in \{+1, -1\}$で、それぞれデータが$C_1, C_2$のどちらに含まれるかを示すとする。係数ベクトルを$w=(w0, w_1, ..., w_d)^T$, 特徴量ベクトルを$x=(x_0=1, x_1, ..., x_d)^T$と表すとすると、

線形識別関数は

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

と表すことができる。

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

../_images/f927c534561a10d7a187e090c62fc275b87b061a92a82c296198006a48a6a666.png

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

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

../_images/26efd72e3f8477609fcbe003898f45f72edb1de909934a8b7b2626750da241ed.png

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

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

訓練データと識別超平面との距離のことをその識別超平面のマージン（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/2ab8c83b591009dd036ebe6147429bdfc94dd29cb147a4cbc8cd08026a479fdf.png

ハードマージンSVM#

学習データの集合を$\mathcal{D}_L = \{(t_i, x_i)\}(i=1,...,N)$とする。係数ベクトルはバイアス項$b$を外に出す形で、$w=(w_1, ..., w_d)^T$と表記する。特徴量ベクトルは$x=(x_1,...,x_d)^T$である。$t_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} t_i = +1\\ w^T x_i + b \leq -1 & \text{if} \hspace{0.5em} t_i = -1 \end{cases} \end{split}\]

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

\[ t_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||}$となる。

最適識別超平面#

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

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

の下で係数ベクトルのノルム$||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} t_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 (t_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 t_i x_i = 0$

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

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

(4) $\alpha_i \geq 0$

(5) $\alpha_i (t_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 t_i w_*^T x_i - b \sum^N_{i=1} \alpha_i t_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 t_i = 0)\\ &= \sum^N_{i=1} \alpha_i - \frac{1}{2} \sum^N_{i=1} \sum^N_{j=1} \alpha_i \alpha_j t_i t_j x_i^T x_j \end{align} \end{split}\]

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

KKT条件(1)より最適解は$w_* = \sum^N_{i=1} \alpha_i t_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 t = 0 \end{align} \end{split}\]

ここで

\[\begin{split} \begin{align} \boldsymbol{1} &= (1,...,1)^T\\ H &= (H_{ij} = t_i t_j x_i^T x_j)\\ t &= (t_1,...,t_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 t \]

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

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

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

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

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

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

実装例（cvxpy）#

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

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

t = y
# Hを作成
n = X.shape[0]
H = np.zeros(shape=(n, n))
for i in range(n):
    for j in range(n):
        H[i, j] = t[i] * t[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 @ t == 0])
prob.solve()
print("The optimal value is", prob.value)
print("alpha is", alpha.value)

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

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

../_images/61c670edb28df5cd14c28fe63f39a46828ee1d3cfd67a57ba1b95e44adf574ab.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/f33a61fdbb9c3536cab52ae757ddcad03bff5add80c03ef3139b2f6000ef618e.png

ソフトマージンSVM#

C-SVM#

# TODO

ν-SVM#

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

カーネルトリック#

# 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 0x7fdc07dc8700>

../_images/61ce6a7dae09b65af86d09b5f7f9cf69bd78dce79f1f180d2c53bd799c933674.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 0x7fdc07bb2290>

../_images/4a8dfd778b21c12c915585855961ba534c7c2a3f17d4d672eea67d877786a6a8.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