サポートベクターマシンの導出からPythonによる2次計画ソルバの利用まで。 サポートベクターマシンとPythonによる2次計画問題：基底展開とカーネル関数の導入に続いている。
まずは2クラス分類から。
2クラス分類
分類（classification）の中でも最も単純なものが2クラス分類になる。\(n\)個の特徴とラベルの組からなるデータセット\(\{(x_i,y_i)\}_{i=1}^n\)から特徴とラベルに関するルールを学習して未知の特徴に対するラベルを正しく予想するのが目的になる。
特徴は\(d\)次元ベクトル、ラベルは\(\{+1,-1\}\)の2種類とし、識別関数\(f(x)\)を持ってきて
\[ y=\mathrm{sign}\;f(x)\]
という形でもってその特徴のラベルを与える識別器（classifier）を考える場合が多い。\(\mathrm{sign}\)は符号を返す関数なので、\(f(x)\ge 0\)であれば\(y=+1\), \(f(x)\lt 0\)であれば\(y=-1\)となる。\(f(x)=0\)の場合をどちらに分類するかという問題は本質的ではない。
\(f(x)=0\)を満たすような\(x\)の集合が決定境界（decision boundary）を与えているので\(f(x)=0\)がどういう集合（直線、平面、曲面など）でどこにあるのかが重要だ。
ハードマージン最大化
\(f(x)\)として次の形式を考える。
\[ f(x)=w^\top x+b\]
ここに\(w\in\mathbb{R}^d\), \(b\in\mathbb{R}\)である。この\(w, b\)をマージン最大化するように選ぶ。マージンというのは決定境界に一番近いデータ点までの距離のことだ。データ点が（超）平面によって完全に分けられると仮定する（線形分離可能という）と
\[ \begin{align} y_i=+1\;\Rightarrow\; w^\top x_i+b\ge0\\ y_i=-1\;\Rightarrow\; w^\top x_i+b\lt0 \end{align}\]
を全ての\(i=1,\dots,n\)で満たすような\(w,b\)がいくらか存在していてその中から\(w^*,b^*\)を選ぶ。それぞれのラベルの中で決定境界に最も近いものを\(x_+^*,x_-^*\)とすると、そのような\(x_+^*,x_-^*\)は
\[ \begin{align} w^\top x_+^*+b&\le w^\top x_i+b &(i=1,\dots,n,\;y_i=+1)\\ w^\top x_-^*+b&\ge w^\top x_i+b &(i=1,\dots,n,\;y_i=-1) \end{align}\]
を満たす。\(w^\top x+b=0\)が境界なので最も際どいものどおしと考えればいい。今考えたいマージンはこのような点だけを用いて計算することができて、
\[ \begin{align} \text{margin} &= \frac{w^\top}{\|w\|}x_+^*-\frac{w^\top}{\|w\|}x_-^*\\ &=\frac{w^\top}{\|w\|}(x_+^*-x_-^*) \end{align}\]
となる。マージンの大きさは点\(x_+^*,x_-^*\)と\(w\)の向きによってのみ決まることがわかる。\(w^\top x+b=0\)が定める決定境界は\(w\)や\(b\)を同じ大きさだけ定数倍しても変わらないし従って\(x_+^*,x_-^*\)の選択にも影響しない。そこで、一般性を損なわずに
\[ \begin{align} w^\top x_+^*+b&=+1\\ w^\top x_-^*+b&=-1 \end{align}\]
という制約を設ける。このとき、
\[ \begin{align} \text{margin}&=\frac{(w^\top x_+^*+b)-(w^\top x_-^*+b)}{\|w\|}\\ &=\frac{2}{\|w\|} \end{align}\]
となり、制約は
\[ \begin{align} +1&\le w^\top x_i+b&\quad(i=1,\dots,n,\; y_i=+1)\\ -1&\ge w^\top x_i+b&\quad(i=1,\dots,n,\; y_i+-1) \end{align}\]
になる。これは\(y_i=\pm 1\)をうまく利用すると
\[ y_i(w^\top x_i+b)\ge 1\quad(i=1,\dots,n)\]
とシンプルに書ける。\(2/\|w\|\)の最大化は\(\|w\|^2\)の最小化と同じことなので、マージン最大化規準による分類は結局
\[ \min_{w,b}\;\frac12\|w\|^2\quad\text{s.t.}\quad y_i(w^\top x_i+b)\ge 1\quad(i=1,\dots,n)\]
の解\(w^*,b^*\)による分類器によるものといえる。最初の図にそのような\(w^*,b^*\)の例を示している。
ソフトマージン最大化
ハードマージン最大化ではデータ点が線形分離可能であることを仮定していた。この場合データ点が線形分離不可能な配置で分布していた場合に最適化問題の解なしということになってしまう。以下の左図などがそうだ。この場合はどんな直線を引いてもオレンジと青を分離することは出来ない。
そこで、スラック変数\(\xi_i\ge0\), \(i=1,\dots,n\)を導入してハードマージンの問題の制約を少しだけなら破ってもいいようにする。
\[ y_i(w^\top x_i+b)\ge 1-\xi_i\quad(i=1,\dots,n)\]
しかしあくまで破っていいのは少しだけ。\(\xi_i\)がいくらでも大きくなってしまうと制約が無いのと同じになってしまうので
\[ \min_{w,b}\;\frac12\|w\|^2+C\sum_{i=1}^m\xi_i\]
というように目的関数を修正する。\(C\)はバランスパラメータだ。以上を踏まえるとソフト―マージンによるマージン最大化問題は次になる。
\[ \begin{align} \min_{w,b}\;\frac12\|w\|^2+C\sum_{i=1}^n&\xi_i\\ \text{s.t.}&\quad \left\{\begin{matrix} y_i(w^\top x_i+b)\ge 1-\xi_i\\ \xi_i\ge 0 \end{matrix}\right.\quad(i=1,\dots,n) \end{align}\]
図にに\(\xi_i\)の意味も示した。\(\xi_i=0\)はマージンの外で余裕を持って分類されている、\(0\lt\xi_i\le1\)はマージンの中にいるが違う側には分類されていない。\(\xi_i\ge1\)は誤った側に分類されている。
ここからはおまけ。条件式の部分は次のように書き換えられる。
\[ \xi_i\ge 1-y_i(w^\top x_i+b)\ge0\quad(i=1,\dots,n)\]
目的関数は\(\xi_i\)を小さくすればするほど良くなる。だから最適値は常に
\[ \begin{align} \xi_i&=\left\{\begin{matrix} 1-y_i(w^\top x_i+b) & 1-y_i(w^\top x_i+b)\ge0\\ 0 & 1-y_i(w^\top x_i+b)\lt0 \end{matrix}\right.\\ &=\max\{0, 1-y_i(w^\top x_i +b)\},\quad(i=1,\dots,n) \end{align}\]
と書ける。これを目的関数に代入して条件を消去すると
\[ \min_{w,b}\;\frac12\|w\|^2+C\sum_{i=1}^n\max\{0, 1-y_i(w^\top x_i+b)\}\]
\(C=1/2\lambda\)としてヒンジ損失\(\ell(y,f)=\max\{0,1-yf\}\)を用いるとこれは
\[ \min_{w,b}\;\sum_{i=1}^n \ell(y_i,w^\top x_i+b)+\lambda\|w\|^2\]
という正則化付きのヒンジ損失最小化によるによる識別関数の決定問題になる。
ヒンジ損失は分類問題における0-1損失を上から抑える形の損失で、原点で不連続な0-1損失の緩和になっている。先の図における\(w^*\)方向に見るとヒンジ損失の値とスラック変数の値が対応していることがわかる。
2次計画法
普通ならここでカーネル化の話に入るが、このままでもPythonの2次計画ソルバで解けることを見る。
まず2次計画法（quadratic programming）は次のような形式に帰着できる最適化問題のこと。
\[ \begin{align} \min_x\;\frac12x^\top Px&+q^\top x\\ &\text{subject to }\left\{\begin{matrix} Gx\le h\\ Ax = b \end{matrix}\right. \end{align}\]
\(x\in\mathbb{R}^d\)のとき、\(P\)は\(d\times d\)行列、\(q\)は\(d\)次元ベクトルになる。不等式制約が\(n\)個あるとき\(G\)は\(n\times d\)行列で\(h\)は\(n\)次元ベクトル、ここでの不等式は各行ごとの比較を表している。\(A,b\)についても等式制約\(m\)個のときそれぞれ\(m\times d\)行列と\(m\)次元ベクトル。
2次計画法の定義を言葉で言うと、目的関数が各変数に対する2次の多項式から成っていて制約は各変数について線形であるもの。2次の多項式は一般的に行列\(P,q\)を用いた上のような2次形式で表せる。\(P\)が正定値ならこの問題の大域的最適解は一意に定まり、計算量的にも効率的に解を求めるアルゴリズムが存在する。
いま、ハードマージンとソフトマージンの最大化問題をこの形式に問題に当てはめてみる。まずハードマージンの場合は簡単で\(x=(w;b)\)、（\(w\)と\(b\)をつなげた\(d+1\)次元ベクトル、説明変数の\(x\)と混同してはいけない！）としてブロック行列を使って、
\[ P=\begin{bmatrix} I_d & \boldsymbol{0}_d\\ \boldsymbol{0}_d^\top & 0 \end{bmatrix},\quad q=\boldsymbol{0}_{d+1},\quad G=-Y(X,\boldsymbol{1}_n),\quad h=-\boldsymbol{1}_n\]
とすればハードマージン最大化の問題を2次計画問題に当てはめられる。等式制約は無い。表記として\(\boldsymbol{0}_d\)や\(\boldsymbol{1}_n\)はそれぞれ\(d\)次元ゼロベクトルと\(n\)次元の要素が全て1のベクトルを表している。少しだけPythonを前提とした書き方になるが\((;)\)は行列を縦につなげる操作（np.r_）、\((,)\)は行列を横に繋げる操作（np.c_）を表している。\(X\)は説明変数\(n\)個を縦に並べた\(n\times d\)行列、\(Y\)はラベル\(\{+1,-1\}\)を縦に並べた\(n\)次元ベクトル。ここでは特別に\(Y(\cdot)\)は行ごとの積を表すことにする。
以上の行列はnumpyを使ってストレートに記述できる。cvxoptを用いてこの2次計画問題を解くコードを示す。
import numpy as np
import cvxopt as opt
# X, Y はそれぞれn,d行列、n,1行列のデータセットとする
n, d = X.shape
P = opt.matrix(np.diag(np.r_[np.ones(d), 0]))
q = opt.matrix(np.zeros(d+1))
G = opt.matrix(-Y.reshape(-1, 1)*np.c_[X, np.ones((n, 1))])
h = opt.matrix(-np.ones(n))
sol = opt.solvers.qp(P, q, G=G, h=h)
x = np.array(sol['x']).ravel()
w, b = x[:d], x[d] #解cvxoptでは独自のmatrixクラスによって行列を表現する。基本的にはnumpyの2次元配列をそのまま渡せばよいが。diagの部分はcvxoptにあるspdiagを用いて次のようにも書ける。
P = opt.spdiag([*np.ones(d), 0])\(P\)は対角成分以外がゼロのスパースな行列なのでこのようにしてcvxoptにスパースであることを伝えた方が内部でよしなにやってくれる可能性があるので場合によってはよい。
次にソフトマージン最大化の場合について。この場合は\(x=(w;b;\xi)\)とするので\(x\)は\(n+d+1\)次元ベクトルになる。
\[ P=\begin{bmatrix} I_d & O_{d\times (n+1)}\\ O_{(n+1)\times d}&O_{(n+1)\times(n+1)} \end{bmatrix}, \quad q=\begin{bmatrix} \boldsymbol{0}_{d+1}\\ C\boldsymbol{1}_{n} \end{bmatrix}\quad G=-\begin{bmatrix} Y(X,\boldsymbol{1}_n) & I_n\\ O_{n\times (d+1)}& I_n \end{bmatrix}, \quad h=-\begin{bmatrix} \boldsymbol{1}_n\\ \boldsymbol{0}_n \end{bmatrix}\]
同じく\(O\)の下付きはそのサイズのゼロ行列だ。何やら複雑になっているがこれでソフトマージン最大化の問題を2次計画問題に当てはめられている。この場合は\(G\)にゼロのブロックやゼロの要素が多いのでcvxoptのスパース表現を用いるのが適切だろう。コード例は以下だ。
import numpy as np
import cvxopt as opt
# X, Y はそれぞれn,d行列、n,1行列のデータセットとする
n, d = X.shape
P = opt.spdiag([*np.ones(d), *np.zeros(1+n)])
q = opt.matrix([*np.zeros(d+1), *(C*np.ones(n))])
G = -opt.sparse([
[
opt.matrix(Y.reshape(-1, 1)*np.c_[X, np.ones((n, 1))]),
opt.spmatrix([], [], [], (n, d+1))
],
[
opt.spdiag([*np.ones(n)]),
opt.spdiag([*np.ones(n)])
]
])
h = -opt.matrix([*np.ones(n), *np.zeros(n)])
sol = opt.solvers.qp(P, q, G=G, h=h)
x = np.array(sol['x']).ravel()
w, b, xi = x[:d], x[d], x[d+1:] #解cvxopt.sparseのブロック行列の並び順はカラムオーダーであることに詰まったので注意。denseとsparseが組み合わさった行列をどのように扱ってくれているのかはわからないが全てdense扱いされるよりは効率を上げる余地があるはず。
実行結果
結果が見やすい2次元の場合において、上に貼ったコードを動かした。データセットを作るところから結果をプロットするところまで含めて以下。
import os, sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as ptcs
import cvxopt as opt
m = 100
rad = -1.2*np.pi/4
A = np.array([
[np.cos(rad), -np.sin(rad)],
[np.sin(rad), np.cos(rad)]
])
# hard-margin classifier
rs = np.random.RandomState(10)
x1 = rs.uniform(-1, 1, m)
y1 = 1/4 + 2*x1**2 + rs.randn(m)*(1+np.abs(x1)*2)/10
X1h = np.c_[x1, y1] @ A.T + 1/4
x2 = rs.uniform(-3/2, 3/2, m)
y2 = -1/4 - x2**2/4 + rs.randn(m)*(1+np.abs(x2)*2)/20
X2h = np.c_[x2, y2] @ A.T + 1/4
X = np.r_[X1h, X2h]
Y = np.r_[np.ones(m), -np.ones(m)]
n, d = X.shape
P = opt.spdiag([*np.ones(d), 0])
q = opt.matrix(np.zeros(d+1))
G = opt.matrix(-Y.reshape(-1, 1)*np.c_[X, np.ones((n, 1))])
h = opt.matrix(-np.ones(n))
sol = opt.solvers.qp(P, q, G=G, h=h)
x = np.array(sol['x']).ravel()
w, b = x[:d], x[d]
xh = np.linspace(-2, 2, 10)
yh = -(w[0]*xh+b)/w[1]
# soft-margin classifier
rs = np.random.RandomState(10)
x1 = rs.uniform(-1, 1, m)
y1 = 1/4 + 2*x1**2 + rs.randn(m)*(1+np.abs(x1)*2)/5
X1s = np.c_[x1, y1] @ A.T + 1/4
x2 = rs.uniform(-3/2, 3/2, m)
y2 = -1/4 - x2**2/4 + rs.randn(m)*(1+np.abs(x2)*2)/10
X2s = np.c_[x2, y2] @ A.T + 1/4
X = np.r_[X1s, X2s]
Y = np.r_[np.ones(m), -np.ones(m)]
C = 10
n, d = X.shape
P = opt.spdiag([*np.ones(d), *np.zeros(1+n)])
q = opt.matrix([*np.zeros(d+1), *(C*np.ones(n))])
G = -opt.sparse([
[
opt.matrix(Y.reshape(-1, 1)*np.c_[X, np.ones((n, 1))]),
opt.spmatrix([], [], [], (n, d+1))
],
[
opt.spdiag([*np.ones(n)]),
opt.spdiag([*np.ones(n)])
]
])
h = -opt.matrix([*np.ones(n), *np.zeros(n)])
sol = opt.solvers.qp(P, q, G=G, h=h)
x = np.array(sol['x']).ravel()
w, b, xi = x[:d], x[d], x[d+1:]
xs = np.linspace(-2, 2, 10)
ys = -(w[0]*xs+b)/w[1]
try:
fig, ax = plt.subplots(1,2,figsize=(8,4.8))
tab10 = plt.get_cmap('tab10')
ax[0].set_title('hard-margin classifier')
ax[0].plot(xh, yh, color=tab10(2), label=r'$w^\top x+b=0$')
ax[0].scatter(X1h[:,0], X1h[:,1], s=3, color=tab10(0), label=r'$y_i=+1$')
ax[0].scatter(X2h[:,0], X2h[:,1], s=3, color=tab10(1), label=r'$y_i=-1$')
ax[0].legend()
ax[1].set_title('soft-margin classifier')
ax[1].plot(xs, ys, color=tab10(2), label=r'$w^\top x+b=0$')
ax[1].scatter(X1s[:,0], X1s[:,1], s=3, color=tab10(0), label=r'$y_i=+1$')
ax[1].scatter(X2s[:,0], X2s[:,1], s=3, color=tab10(1), label=r'$y_i=-1$')
ax[1].legend()
fig.tight_layout()
plt.show()
finally:
plt.close('all')
pcost dcost gap pres dres
0: 2.5414e-01 2.1188e+02 9e+02 2e+00 3e+02
1: 2.9352e+00 -1.0316e+02 4e+02 9e-01 1e+02
2: 5.1157e+00 -1.2127e+02 4e+02 8e-01 1e+02
3: 8.6748e+00 -5.1909e+01 4e+02 7e-01 9e+01
4: 1.9873e+01 -1.3486e+02 4e+02 5e-01 7e+01
5: 4.6858e+01 -6.4396e+01 1e+02 8e-02 1e+01
6: 4.8872e+01 -5.6007e+01 1e+02 6e-02 9e+00
7: 4.5801e+01 1.8120e+01 3e+01 2e-15 5e-14
8: 4.4532e+01 3.4773e+01 1e+01 2e-15 4e-14
9: 4.3396e+01 4.3116e+01 3e-01 2e-15 1e-13
10: 4.3352e+01 4.3349e+01 3e-03 3e-15 1e-13
11: 4.3352e+01 4.3352e+01 3e-05 2e-15 1e-13
Optimal solution found.
pcost dcost gap pres dres
0: -1.5361e+04 1.6432e+04 6e+04 3e+01 8e+00
1: 3.3254e+03 -3.4728e+03 1e+04 4e+00 1e+00
2: 8.9616e+02 -4.6765e+02 2e+03 5e-01 2e-01
3: 4.3336e+02 8.0559e+01 5e+02 8e-02 2e-02
4: 3.3711e+02 1.5291e+02 2e+02 3e-02 1e-02
5: 3.1074e+02 1.7879e+02 2e+02 2e-02 7e-03
6: 2.7781e+02 2.0605e+02 8e+01 8e-03 2e-03
7: 2.6940e+02 2.1422e+02 6e+01 6e-03 2e-03
8: 2.5613e+02 2.2556e+02 3e+01 2e-03 7e-04
9: 2.4707e+02 2.3294e+02 2e+01 8e-04 3e-04
10: 2.4226e+02 2.3739e+02 5e+00 8e-05 2e-05
11: 2.3997e+02 2.3927e+02 7e-01 1e-05 3e-06
12: 2.3961e+02 2.3958e+02 2e-02 2e-07 6e-08
13: 2.3959e+02 2.3959e+02 2e-04 2e-09 6e-10
Optimal solution found.まとめ
2次計画問題のソルバがscipyにありそうでない。cvxoptを用いてサポートベクターマシンの最適化を素朴に解いた。 サポートベクターマシンとPythonによる2次計画問題：基底展開とカーネル関数の導入に続く。
参考文献
福水健次著「カーネル法入門―正定値カーネルによるデータ解析―」、朝倉書店2016. pp. 47, 3.4「サポートベクターマシンの基本事項」
CVXOPT User's Guide (Quadratic programming) https://cvxopt.org/userguide/coneprog.html#quadratic-programming