Deploy on Friday

線形SVMの理論と実装

December 06, 2018

このエントリは機械学習の数理 Advent Calendar 2018の7日目の記事になります。ハードマージン線形SVMを一部ソルバのちからを借りつつPythonで実装します。

数学的な準備

超平面について

pp次元の分離超平面は、

w1x1+w2x2+⋯+wpxp+b=0w_1x_1 + w_2x_2 + \dots + w_px_p + b = 0

という式で表されます。 ベクトル

w=(w1,w2,…,wp)Tx=(x1,x2,…,xp)T\begin{aligned} \boldsymbol{w} &= (w_1, w_2, \dots, w_p)^{\mathsf{T}}\\ \boldsymbol{x} &= (x_1, x_2, \dots, x_p)^{\mathsf{T}} \end{aligned}

を用いて

wx+b=0\boldsymbol{wx} + b = 0

と表されます。 また、超平面を実数倍しても、超平面自体は変わりません。

rwx+rb=0r\boldsymbol{wx} + rb = 0

超平面と点の距離

高校生の時に、1次関数(w1x1+w2x2+b=0w_1x_1 + w_2x_2 + b = 0)と点(x1′,x2′)(x_{1}', x_{2}')の距離の公式を

d=∣w1x1′+w2x2′+b∣w12+w22d = \frac{|w_1x_1' + w_2x_2' + b|}{\sqrt{w_1^2 + w_2^2}}

と習ったと思います。

dist

これを一般化すると、

d=∣w1x1′+w2x2′+,…,+wpxp′+b∣w12+w22+⋯+wp2=∣wTx′+b∣∣∣w∣∣\begin{aligned} d &= \frac{|w_1x_1' + w_2x_2' +, \dots, + w_px_p' + b|}{\sqrt{w_1^2 + w_2^2 + \dots + w_p^2}}\\ & = \frac{|\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}' + b|}{||\boldsymbol{w}||} \end{aligned}

となります。

ラグランジュ理論

不等式制約つき最適化問題

min⁡w,bf(w)subject to gi(w)≤0 (i=1,…,n)\min_{\boldsymbol{w}, b} f(\boldsymbol{w})\\ \text{subject to}\ g_i(\boldsymbol{w}) \leq 0\ (i=1, \dots, n)

について

L(w,α)=f(w)+∑i=1nαigi(w)L(\boldsymbol{w, \alpha}) = f(\boldsymbol{w}) + \sum_{i=1}^n \alpha_ig_i(\boldsymbol{w})

をラグランジュ関数と定めたとき、次の4つの条件を満たすw\boldsymbol{w}とα\boldsymbol{\alpha}を見つける問題に帰着される。

(1)∂∂wjL(w,α)=∂∂wjf(w)+∑i=1nαj∂∂wjgi(w)=0(2)∂∂αiL(w,α)=gi(w)≤0(3)αi≥0(4)αigi(w)=0\begin{aligned} (1) & \frac{\partial}{\partial w_j}L(\boldsymbol{w, \alpha}) = \frac{\partial}{\partial w_j}f(\boldsymbol{w}) + \sum_{i=1}^n \alpha_j \frac{\partial}{\partial w_j}g_i(\boldsymbol{w}) = 0\\ (2) & \frac{\partial}{\partial \alpha_i} L(\boldsymbol{w, \alpha})=g_i(\boldsymbol{w}) \leq 0\\ (3) & \alpha_i \geq 0\\ (4) & \alpha_ig_i(\boldsymbol{w}) = 0 \end{aligned}

これらの条件のことをKKT条件(カルーシュ・クーン・タッカー条件)という。詳細については略。

ハードマージン線形SVM

2次計画問題への落し込み

線形分離可能とは

linearly separable

とした時に、

yi(wTxi+b)>0  (i=1,…,n)y_i( \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b) > 0\ \ (i=1,\dots,n)

を満たすことを言う。

support vector

この時、超平面から最も近いデータをサポートベクターと呼び、そのサポートベクターと超平面の距離ddを最大化することを考えます。(ここで超平面と点の距離の公式を使います。)

max⁡w,bd    subject to  yi(wTxi+b)∣∣w∣∣≥d  (i=1,…,n)\begin{aligned} \max_{\boldsymbol{w}, b}d\ \ \ \ \text{subject to}\ \ \frac{y_i(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b)}{||\boldsymbol{w}||} \geq d \ \ (i=1,\dots,n) \end{aligned}

しかし、最大化問題は難しい… そこで、先ほどの制約条件式の両辺をddで割って、

yi(1d∣∣w∣∣wTxi+1d∣∣w∣∣b)≥1  (i=1,…,n)\begin{aligned} y_i\Biggl(\frac{1}{d||\boldsymbol{w}||}\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + \frac{1}{d||\boldsymbol{w}||}b\Biggr) \geq 1 \ \ (i=1,\dots,n) \end{aligned}

超平面を実数倍しても超平面自体は変わらないので、スケールを変換して、

yi(w′Txi′+b)≥1  (i=1,…,n)y_i( \boldsymbol{w}'^{T}\boldsymbol{x}'_i + b) \geq 1\ \ (i=1,\dots,n)

という式が出てきますが、こちらが制約条件式になります。 この時、サポートベクターについては次式が成り立つ。

yi(w′Txi′+b)=1y_i( \boldsymbol{w}'^{T}\boldsymbol{x}'_i + b) = 1

yiy_iは11or−1-1なので、サポートベクターと超平面の距離d′d'は

d′=∣w′Tx′+b∣∣∣w′∣∣=1∣∣w′∣∣\begin{aligned} d' = \frac{|\boldsymbol{w}'^{\mathsf{T}}\boldsymbol{x}' + b|}{||\boldsymbol{w}'||} = \frac{1}{||\boldsymbol{w}'||} \end{aligned}

となります。

dist2

以降は′'を外して説明します。 距離ddを最大化するという問題を解きやすくするために逆数の2乗の最小化問題に置き換えます。 こうすることによって2次計画問題(解が一意に決定する。)に落としこむことができました。

min⁡w,b12∣∣w∣∣2    subject to   yi(wTxi+b)≥1  (i=1,…,n)\begin{aligned} \min_{\boldsymbol{w},b} \frac{1}{2}||\boldsymbol{w}||^2\ \ \ \ \text{subject to}\ \ \ y_i(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b) \geq 1\ \ (i=1,\dots,n) \end{aligned}

主問題から双対問題への落し込み

この不等式制約つき最適化問題をラグランジュ理論を基に双対問題へ落し込む。制約条件式を変形して、

yi(wTxi+b)≥1  (i=1,…,n)−{yi(wTxi+b)−1}≤0  (i=1,…,n)\begin{aligned} y_i(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b) &\geq 1\ \ (i=1,\dots,n)\\ -\{y_i(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b) - 1\}&\leq 0\ \ (i=1,\dots,n) \end{aligned}

ラグランジュ未定乗数 α=(α1,α2,…,αn)T  (αi≥0)\boldsymbol{\alpha}=(\alpha_1, \alpha_2, \dots, \alpha_n)^{\mathsf{T}} \ \ (\alpha_i \geq 0)を定義すると、ラグランジュ関数は、

L(w,b,α)=12∣∣w∣∣2−∑i=1nαi{yi(wTxi+b)−1}\begin{aligned} L(\boldsymbol{w}, b, \boldsymbol{\alpha}) = \frac{1}{2}||\boldsymbol{w}||^2 - \sum_{i=1}^{n}\alpha_i\{y_i(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b) - 1\} \end{aligned}

と表され、偏微分すると、

∂∂bL(w,b,α)=−∑i=1nαiyi=0↔∑i=1nαiyi=0\begin{aligned} \frac{\partial}{\partial b}L(\boldsymbol{w}, b, \boldsymbol{\alpha}) = - \sum_{i=1}^{n}\alpha_iy_i = 0 \leftrightarrow \sum_{i=1}^{n}\alpha_iy_i = 0 \end{aligned} ∂∂wL(w,b,α)=w−∑i=1nαiyixi=0↔w=∑i=1nαiyixi\begin{aligned} \frac{\partial}{\partial \boldsymbol{w}}L(\boldsymbol{w}, b, \boldsymbol{\alpha}) = \boldsymbol{w} - \sum_{i=1}^{n}\alpha_iy_ix_i = 0 \leftrightarrow \boldsymbol{w} = \sum_{i=1}^{n}\alpha_iy_ix_i \end{aligned}

∑i=1nαiyi=0\sum_{i=1}^{n}\alpha_iy_i = 0とw=∑i=1nαiyixi\boldsymbol{w} = \sum_{i=1}^{n}\alpha_iy_ix_iをL(w,b,α)L(\boldsymbol{w}, b, \boldsymbol{\alpha})に代入して計算すると、

LD(α)=12∣∣w∣∣2−∑i=1nαi{yi(wTxi+b)−1}=12∑i=1nαiαjyiyjxiTxj−∑i=1n{αiyi(wTxi+b)−αi}=12∑i=1nαiαjyiyjxiTxj−∑i=1nαiyiwTxi−b∑i=1nαiyi+∑i=1nαi=∑i=1nαi−12∑i=1n∑j=1nαiαjyiyjxiTxj  (∵∑i=1nαiyi=0)\begin{aligned} L_D(\boldsymbol{\alpha}) & = \frac{1}{2}||\boldsymbol{w}||^2 - \sum_{i=1}^n \alpha_i \{ y_i(\boldsymbol{w}^\mathsf{T}\boldsymbol{x}_i+ b) - 1 \} \\ & = \frac{1}{2}\sum_{i=1}^n\alpha_i\alpha_jy_iy_j\boldsymbol{x}_i^\mathsf{T}\boldsymbol{x}_j - \sum_{i=1}^n\{\alpha_iy_i(\boldsymbol{w}^\mathsf{T}\boldsymbol{x}_i+b)-\alpha_i\}\\ & = \frac{1}{2}\sum_{i=1}^n\alpha_i\alpha_jy_iy_j\boldsymbol{x}_i^\mathsf{T}\boldsymbol{x}_j - \sum_{i=1}^n\alpha_iy_i\boldsymbol{w}^\mathsf{T}\boldsymbol{x}_i -b\sum_{i=1}^n\alpha_iy_i + \sum_{i=1}^n\alpha_i\\ & = \sum_{i=1}^{n}\alpha_i - \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_jy_iy_j\boldsymbol{x}^{\mathsf{T}}_i\boldsymbol{x}_j \ \ (\because \sum_{i=1}^n\alpha_iy_i = 0) \end{aligned}

となる。 従って得られた双対問題は、

max⁡αLD(α)    subject to   ∑i=1nαiyi=0   αi≥0 (i=1,…,n)\begin{aligned} \max_{\boldsymbol{\alpha}} L_D(\boldsymbol{\alpha})\ \ \ \ \text{subject to}\ \ \ \sum_{i=1}^{n}\alpha_iy_i = 0 \ \ \ \alpha_i \geq 0 \ (i=1,\dots,n) \end{aligned}

これを2次計画問題のソルバーを使って解を得ることを考える。Pythonにはcvxoptという数理最適化用ライブラリがあるのでこれの中にあるcvxopt.solvers.qpを用いる。先程得られた双対問題をこれに合わせて変形することを考える。

一般的な凸２次計画問題の最適条件の定式は

min⁡α 12αTQα+pα\min_{\boldsymbol{\alpha}} \ \frac{1}{2} \boldsymbol{\alpha}^{\mathsf{T}} Q \boldsymbol{\alpha} + \boldsymbol{p\alpha} subject to Gα≤h Aα=b\begin{aligned} \text{subject to} &\ G\boldsymbol{\alpha} \leq \boldsymbol{h} \\ &\ A\boldsymbol{\alpha} = \boldsymbol{b} \end{aligned}



と表される。一方、最適化したい問題は

min⁡α12∑i=1n∑j=1nαiαjyiyjxiTxj−∑i=1nαi\min_{\boldsymbol{\alpha}} \frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} \alpha_i\alpha_jy_iy_j\boldsymbol{x}_i^\mathsf{T}\boldsymbol{x}_j - \sum_{i=1}^{n}\alpha_i subject to  αi≥0 (i=1,…,n)  ∑i=1nαiyi=0\begin{aligned} \text{subject to} &\ \ \alpha_i \geq 0 \ (i=1,\dots,n) \\ &\ \ \sum_{i=1}^{n}\alpha_iy_i = 0 \end{aligned}

なのでこれを式変形して

min⁡α12(α1,α2,⋯ ,αn)(y1y1x1Tx1y1y2x1Tx2⋯y1ynx1Txny2y1x2Tx1y2y2x2Tx2⋯y2ynx2Txn⋮⋮⋱⋮yny1xnTx1yny2xnTx2⋯ynynxnTxn)(α1α2⋮αn)+(−1,−1,⋯ ,−1)(α1α2⋮αn)=min⁡α12αT(y1y1x1Tx1y1y2x1Tx2⋯y1ynx1Txny2y1x2Tx1y2y2x2Tx2⋯y2ynx2Txn⋮⋮⋱⋮yny1xnTx1yny2xnTx2⋯ynynxnTxn)⏟Qα+(−1)T⏟pα\min_{\boldsymbol{\alpha}} \frac{1}{2}(\alpha_1, \alpha_2, \cdots, \alpha_n) \left( \begin{array}{ccc} y_1y_1\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_1 & y_1y_2\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_1y_n\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_n \\ y_2y_1\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_1 & y_2y_2\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_2y_n\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_n \\ \vdots & \vdots & \ddots & \vdots \\ y_ny_1\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_1 & y_ny_2\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_ny_n\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_n \end{array} \right)\left( \begin{array}{ccc} \alpha_1\\ \alpha_2\\ \vdots\\ \alpha_n \end{array} \right) + (-1, -1, \cdots, -1)\left( \begin{array}{ccc} \alpha_1\\ \alpha_2\\ \vdots\\ \alpha_n \end{array} \right)\\ = \min_{\boldsymbol{\alpha}} \frac{1}{2}\boldsymbol{\alpha}^{\mathsf{T}} \underbrace{ \left( \begin{array}{ccc} y_1y_1\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_1 & y_1y_2\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_1y_n\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_n \\ y_2y_1\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_1 & y_2y_2\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_2y_n\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_n \\ \vdots & \vdots & \ddots & \vdots \\ y_ny_1\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_1 & y_ny_2\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_ny_n\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_n \end{array} \right) }_{Q} \boldsymbol{\alpha} + \underbrace{(\boldsymbol{-1})^{\mathsf{T}}}_{\boldsymbol{p}}\boldsymbol{\alpha} \\ subject to  αi≥0 (i=1,…,n)\text{subject to} \ \ \alpha_i \geq 0 \ (i=1,\dots,n) \\ (α1α2⋮αn)≥(00⋮0)\left( \begin{array}{ccc} \alpha_1\\ \alpha_2\\ \vdots\\ \alpha_n \end{array} \right) \geq \left( \begin{array}{ccc} 0\\ 0\\ \vdots\\ 0 \end{array} \right) (10⋯001⋯0⋮⋮⋱000⋯1)(α1α2⋮αn)≥(00⋮0)∴−(10⋯001⋯0⋮⋮⋱000⋯1)⏟Gα≤(00⋮0)⏟h\begin{aligned} \left( \begin{array}{ccc} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & \cdots & 1 \end{array} \right) \left( \begin{array}{ccc} \alpha_1\\ \alpha_2\\ \vdots\\ \alpha_n \end{array} \right) & \geq \left( \begin{array}{ccc} 0\\ 0\\ \vdots\\ 0 \end{array} \right) \\ \therefore \underbrace{ - \left( \begin{array}{ccc} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & \cdots & 1 \end{array} \right) }_{G} \boldsymbol{\alpha} & \leq \underbrace{ \left( \begin{array}{ccc} 0\\ 0\\ \vdots\\ 0 \end{array} \right)}_{\boldsymbol{h}} \end{aligned}

また、もう一つの条件式についても

∑i=1nαiyi=0\sum_{i=1}^{n} \alpha_iy_i = 0 (y1,y2,⋯ ,yn)(α1α2⋮αn)=0(y_1, y_2, \cdots, y_n) \left( \begin{array}{ccc} \alpha_1\\ \alpha_2\\ \vdots\\ \alpha_n \end{array} \right) = 0 (y1,y2,⋯ ,yn)⏟Aα=0⏟b\underbrace{(y_1, y_2, \cdots, y_n)}_{A} \boldsymbol{\alpha} = \underbrace{0}_{b}

つまるところ、

Q=(y1y1x1Tx1y1y2x1Tx2⋯y1ynx1Txny2y1x2Tx1y2y2x2Tx2⋯y2ynx2Txn⋮⋮⋱⋮yny1xnTx1yny2xnTx2⋯ynynxnTxn),Q = \left( \begin{array}{ccc} y_1y_1\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_1 & y_1y_2\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_1y_n\boldsymbol{x}_1^\mathsf{T}\boldsymbol{x}_n \\ y_2y_1\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_1 & y_2y_2\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_2y_n\boldsymbol{x}_2^\mathsf{T}\boldsymbol{x}_n \\ \vdots & \vdots & \ddots & \vdots \\ y_ny_1\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_1 & y_ny_2\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_2 & \cdots & y_ny_n\boldsymbol{x}_n^\mathsf{T}\boldsymbol{x}_n \end{array} \right), p=(−1)T,\boldsymbol{p} = (\boldsymbol{-1})^{\mathsf{T}}, G=−(10⋯001⋯0⋮⋮⋱000⋯1),G = - \left( \begin{array}{ccc} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & \cdots & 1 \end{array} \right), h=(00⋮0),\boldsymbol{h} = \left( \begin{array}{ccc} 0\\ 0\\ \vdots\\ 0 \end{array} \right),\\ A=(y1,y2,⋯ ,yn),A = (y_1, y_2, \cdots, y_n),\\ b=0.b = 0.

として2次計画問題を解けば良いことになる。ただし、このときAAは11行nn列の行列として扱う。

これにより得られた解α=(α1,α2,⋯ ,αn)\boldsymbol{\alpha} = (\alpha_1, \alpha_2, \cdots, \alpha_n)とすると、4番目のKKT条件から

αigi(w)=0\alpha_ig_i(\boldsymbol{w}) = 0 αi(yi(wTxi+b)−1)=0\alpha_i \left(y_i(\boldsymbol{w}^\mathsf{T}\boldsymbol{x}_i + b) - 1\right) = 0

を満たすことになるが、サポートベクターについてはwTxi+b=1 or −1\boldsymbol{w}^\mathsf{T}\boldsymbol{x}_i + b = 1\ \text{or}\ -1となるため、αi\alpha_iがどのような値であっても条件を満たすが、それ以外についてはyi(wTxi+b)−1y_i(\boldsymbol{w}^\mathsf{T}\boldsymbol{x}_i + b)-1が00にならないため、必然的にαi\alpha_iが0になる必要がある。このため、w\boldsymbol{w}を求めるには非ゼロのαi\alpha_iについてのみ考えれば良く、

w=∑i∈Sαiyixiwhere S is the set of indices of support vectors.\boldsymbol{w} = \sum_{i \in S}\alpha_iy_i\boldsymbol{x}_i\\ \text{where\ }S \text{ is the set of indices of support vectors.}

と求めることができる。一方バイアスbbについて考える。サポートベクターについては、

yi(wTxi+b)=1y_i\left(\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_i + b\right) = 1

が成り立つので、正例のサポートベクターをx+x_+、負例のサポートベクターをx−x_-とすると、

wTx++b=1wTx−+b=−1↔wTx++wTx−+2b=0↔b=−12(wTx++wTx−)\boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_+ + b = 1\\ \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_- + b = -1\\ \leftrightarrow \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_+ + \boldsymbol{w}^{\mathsf{T}}\boldsymbol{x}_- + 2b = 0\\ \leftrightarrow b = -\frac{1}{2}(\boldsymbol{w^\mathsf{T}x_+} + \boldsymbol{w^\mathsf{T}x_-})

Pythonで実装すると

import numpy as np
import cvxopt
import matplotlib.pyplot as plt
import sklearn.datasets
from typing import List

N = 100
positive: List[np.ndarray] = []
negative: List[np.ndarray] = []

mean1 = [-2, 2]
mean2 = [2, -2]
cov = [[1.0, 0.3], [0.3, 1.0]]
    
positive.extend(np.random.multivariate_normal(mean1, cov, N//2))
negative.extend(np.random.multivariate_normal(mean2, cov, N//2))
    
X = np.vstack((positive, negative))
    
y = []
for i in range(N//2):
    y.append(1.0)   # 正例 y_i = 1
for i in range(N//2):
    y.append(-1.0)  # 負例 y_i = -1
y = np.array(y)


_Q = np.zeros((N, N))
for i in range(N):
    for j in range(N):
        _Q[i, j] = y[i]*y[j]*np.dot(X[i], X[j])
Q = cvxopt.matrix(_Q)
p = cvxopt.matrix(-np.ones(N))
G = cvxopt.matrix(-np.eye(N))
h = cvxopt.matrix(np.zeros(N))
A = cvxopt.matrix(y[np.newaxis], (1, N), 'd')
b = cvxopt.matrix(0.0)
solution = cvxopt.solvers.qp(Q, p, G, h, A, b)

alpha = np.array(solution['x']).flatten()

top2_sv_indices = alpha.argsort()[-2:]
sv_indices = alpha > 1e-6
W = np.dot(alpha[sv_indices] * y[sv_indices], X[sv_indices])
bias = - np.dot(X[top2_sv_indices], W).mean()

xs = np.array([X.min(), X.max()])
ys = np.array([(-W[0]*xs[0]-bias)/W[1], (-W[0]*xs[1]-bias)/W[1]])


plt.scatter(X[:, 0][y == 1], X[:, 1][y == 1], label="positive")
plt.scatter(X[:, 0][y == -1], X[:, 1][y == -1], label="negative")
for coordinate in X[sv_indices]:
    plt.annotate('sv', coordinate)
plt.plot(xs, ys)
plt.legend()
plt.grid()
plt.show()

実行した出力

     pcost       dcost       gap    pres   dres
 0: -7.4843e+00 -1.2390e+01  3e+02  2e+01  2e+00
 1: -4.4900e+00 -1.3805e+00  3e+01  2e+00  2e-01
 2: -1.0676e-01 -4.8913e-01  8e-01  2e-02  2e-03
 3: -1.8470e-01 -2.8934e-01  1e-01  3e-03  3e-04
 4: -2.4603e-01 -2.7739e-01  3e-02  6e-05  6e-06
 5: -2.6911e-01 -2.7311e-01  4e-03  6e-06  7e-07
 6: -2.7247e-01 -2.7268e-01  2e-04  2e-07  2e-08
 7: -2.7266e-01 -2.7266e-01  2e-06  2e-09  2e-10
 8: -2.7266e-01 -2.7266e-01  2e-08  2e-11  2e-12
 Optimal solution found.
hard svm

おわりに

今回は完全に線形分離可能な問題下におけるハードマージン線形SVMについて扱いました。線形分離可能でない場合については、今回解いた2次計画問題の設定では解にたどり着けないのでスラック変数というのを導入するのですが、それに関してはまた次回以降。

参考

Written by@Minato Sato
Senior Software Engineer - Embedded AI

GitHubTwitterFacebookLinkedIn

© 2023 Minato Sato. All Rights Reserved