サポートベクターマシンとPythonによる2次計画問題：導入からスラック変数までからの続き。
サポートベクターマシン=カーネル法というレベルでサポートベクターマシンにはカーネル法の概念が必須。前段階としての基底展開からカーネル法の2次計画問題のPythonによる計算など。
線形な決定境界：linear boundary
ハードマージンの場合に解なしとならないようにスラック変数導入してソフトマージンにしたはいいが、説明変数に対して線形な決定境界しか生成しないことに変わりない。
こんな感じで線形分離不可能なデータ分布に対してはいくら\(C\)というパラメータを調整しても有効な境界を与えそうにない。データの分布は見た目からしてこんなに離れているにもかかわらずだ。このデータは2次関数とかを用いればいくらでも簡単に分離出来るはずなのに。
（もっとも、この例だとこういう分け方もすごく悪いわけではないんだがそれは置いとく。）
基底展開：basis expansion
そこでカーネル法という無限個の基底展開に相当するものでサポートベクターマシンを作り直すところだが一旦基底展開についてみる。
2次元の場合にしか手にを得ないので2次元の場合を考える。このとき\(x=(x_1,x_2)\)という特徴を関数\(\Phi:\mathbb{R}^2\to\mathbb{R}^d\)によって\(d\)次元の特徴に展開する。ちなみに、一般に\(\Phi:\mathbb{R}^p\to\mathbb{R}^d\)という場合に
- \(p\gt d\)のとき特徴抽出
- \(p=d\)のとき特徴変換
- \(p\lt d\)のとき基底展開
などという。線形写像が同次元のときに線形変換というような感じ。最後は特徴展開でもいいような気はするがそれはいい。今考えるのは
\[ \Phi(x)=\Phi(x_1,x_2):=\begin{bmatrix} x_1\\ x_2\\ x_1^2\\ x_1x_2\\ x_2^2 \end{bmatrix}\in\mathbb{R}^5\]
だ。ようは2次の項をもともとの特徴に追加した。今までの\(x\)を\(\Phi(x)\)に置き換えて同じようにソフトマージンで上に図のデータセットを分類してみる。するとこう。
このように曲がった決定境界を表現することが出来る。この場合は全てのデータを完全に分離できているので\(\xi_i=0, i=1,\dots,n\)でハードマージンの場合に一致する。決定境界が非線形となるためには、非線形な特徴で基底展開をしているかどうかが重要だ。少し違う例でやってみよう。
\[ \Phi(x)=\Phi(x_1,x_2):=\begin{bmatrix}x_1\\x_2\\\cos(x_1x_2)\end{bmatrix}\]
とするとこう。
このように、2次の項全部を追加しなくても何かしら非線形な特徴を用いれば非線形な決定境界を定めることができる。ここでは2次元の説明変数の場合に注目したが、例えば\(d\)次元の場合にこれをやるとなるとどうなるだろうか。単に非線形な特徴を展開するといっても\(d\)個も変数があれば、ここで見た\(\cos\)を入れるなども視野に入れ始めると選び方はほぼ無限大になってしまう。
ちなみに基底展開という話はカーネル法がそうであるようにサポートベクターマシンに限った話ではない。
カーネル関数
そこでカーネル法。サポートベクターマシン=カーネル法の代表といって違いない。基底展開を含むサポートベクターマシンの中で内積を使っていた\(w^\top\Phi(x)\)の部分に注目する。カーネル法の本質のひとつにはこのような内積の演算を無限次元に拡張することがある。
基底を無限個に展開していくという操作を考えると、
\[ w=\begin{bmatrix}w_1\\w_2\\\vdots\end{bmatrix},\quad\Phi(x)=\begin{bmatrix}\Phi_1(x)\\\Phi_2(x)\\\vdots\end{bmatrix}\]
というように\(w\)と\(\Phi\)の値域は無限に縦長なベクトルになる。無限に縦長なベクトルと関数のアナロジーから、この2つの内積演算をRKHS（再生核ヒルベルト空間）上の内積演算に置き換えることは自然に考えれる。
\[ \begin{gather}w\to h(\cdot),\;\Phi:\mathbb{R}^d\to\mathcal{H},\\\Phi:x\mapsto k_x(\cdot)\\w^\top\Phi(x)\to\langle h, k_x\rangle_\mathcal{H}\end{gather}\]
という対応になる。同時に\(\|w\|^2=w^\top w\to\langle h, h\rangle_\mathcal{H}=\|h\|_\mathcal{H}\)も考える。これを踏まえてソフトマージンの最大化問題を書き直すと次になる。
\[ \begin{align}\min_{h\in\mathcal{H},b,\xi}\;\frac12&\|h\|_\mathcal{H}+C\sum_{i=1}^n\xi_i\\&\text{s.t. }\left\{\begin{matrix}y_i(\langle h,k_{x_i}\rangle_\mathcal{H}+b)\ge 1-\xi_i\\\xi_i\ge 0\end{matrix}\right.\quad (i=1,\dots,n)\end{align}\]
対応通りに書き直しただけだ。ここでRepresenter定理を使って\(\mathcal{H}\)という関数空間から最適な関数\(h\)を選ぶという、一見してやりようのない問題を再び有限次元ベクトルの最適化に書き直す。
Representer定理のいうところによると、こういう問題の場合考える\(h\)は次の形式だけで充分である。
\[ h(x)=\sum_{j=1}^n\alpha_jk_{x_j}(x)\]
これをもとの問題に代入するのだが、そこで\(\langle k_x,k_{x}'\rangle_\mathcal{H}=k(x,x')\)という性質を用いる。\(k_x\)は\(x\)を無限個に特徴展開した無限次元ベクトルのアナロジーから来る関数で、別な\(x'\)を無限個に特徴展開したのに対応する\(k_{x'}\)との内積計算には一見して
\[ \langle k_x,k_{x'}\rangle_\mathcal{H}=\int k_x(\tau)k_{x'}(\tau)d\tau\]
というような量の近似計算が必要になるとも思えそうなところが、RKHSにはある関数\(k(x,x')\)が存在してどんな\(x,x'\)に対しても\(\langle k_x,k_{x}'\rangle_\mathcal{H}=k(x,x')\)と出来ることが知られている。逆もいえて、\(k(x,x')\)を一つ定めればそれに対応するRKHSが一つ定まる。だから以降は\(k(x,x')\)という関数ありきで進めてRKHSどうこうは考えない（こういうやり方をカーネルトリックともいう）。
\[ \begin{align}\|h\|_\mathcal{H}&=\left\langle\sum_{i=1}^n\alpha_ik_{x_i},\sum_{j=1}^n\alpha_jk_{x_j}\right\rangle_\mathcal{H}\\&=\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_j\langle k_{x_i},k_{x_j}\rangle_\mathcal{H}\\&=\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_jk(x_i,x_j)\\\langle h,k_{x_i}\rangle_\mathcal{H}&=\left\langle\sum_{j=1}^n\alpha_jk_{x_j},k_{x_i}\right\rangle_\mathcal{H}\\&=\sum_{j=1}^n\alpha_j\langle k_{x_j},k_{x_i}\rangle_\mathcal{H}\\&=\sum_{j=1}^n\alpha_jk(x_j,x_i)\end{align}\]
より
\[ \begin{align}\min_{\alpha,b,\xi}\;\frac12\sum_{i=1}^n\sum_{j=1}^n&\alpha_i\alpha_jk(x_i,x_j)+C\sum_{i=1}^n\xi_i\\&\text{s.t.}\left\{\begin{matrix}y_i\left(\sum_{j=1}^n\alpha_jk(x_j,x_i)+b\right)\ge 1-\xi_i\\\xi_i\ge 0\end{matrix}\right. \quad(i=1,\dots,n)\end{align}\]
が考えるべき問題。グラム行列\(K_{ij}=k(x_i,x_j)\)とベクトル\(\alpha=(\alpha_1,\dots,\alpha_n)\)を用いると
\[ \begin{align}\min_{\alpha,b,\xi}\;\frac12\alpha^\top K\alpha+C\sum_{i=1}^n\xi_i\\&\text{s.t.}\left\{\begin{matrix}y_i\left(K\alpha+b\right)\ge 1-\xi_i\\\xi_i\ge 0\end{matrix}\right. \quad(i=1,\dots,n)\end{align}\]
と書ける。これはもともとの説明変数の次元\(d\)によらない\(2n+1\)次元ベクトルの最適化になっている。無限次元の基底展開を扱う計算を高々\(d\)を\(n\)に増やしただけで実現できるなら安いもんだろう？
最終的にこの問題の解\(\alpha^*\)を用いて決定境界は次で与えられる。
\[ f(x)=\sum_{j=1}^n\alpha_j^*k(x_j,x)+b=0\]
これは、\(\boldsymbol{k}(x)=(k(x_1,x),\dots,k(x_n,x))\)というベクトルを用いると
\[ f(x)=\boldsymbol{k}(x)^\top\alpha^*+b=0\]
と簡単に書ける。
カーネル化SVMの2次計画問題
カーネル関数を用いたサポートベクターマシンの最適化問題はほとんど以下の2次計画問題と同じ形をしている。
\[ \begin{align}\min_x\;\frac12x^\top Px&+q^\top x\\&\text{s.t. }\left\{\begin{matrix}Gx\le h\\Ax = b\end{matrix}\right.\end{align}\]
最適化変数を\(x=(\alpha;b;\xi)\)としてブロック行列をつかうと
\[ \begin{gather}P=\begin{bmatrix}K & O_{n\times(n+1)}\\O_{(n+1)\times n} & O_{(n+1)\times(n+1)}\end{bmatrix},\quad q=\begin{bmatrix}\boldsymbol{0}_{(n+1)}\\C\boldsymbol{1}_n\end{bmatrix}\\G=-\begin{bmatrix}Y(K,\boldsymbol{1}_n) & I_n\\O_{n\times(n+1)} & I_n\end{bmatrix},\quad h=-\begin{bmatrix}\boldsymbol{1}_{n}\\\boldsymbol{0}_{n}\end{bmatrix}\end{gather}\]
となる。カーネル関数の正定値性よりグラム行列は正定値なのでこの2次計画問題は大域的最適解が求まる。見てわかるようにこれはソフトマージンの最適化問題での\(X\)を\(K\)に変えただけ。そのままグラム行列に置き換えるだけで動く。実際に計算するときには数値的安定性のために小さい\(\delta\)をつかって\(K\)の代わりに\(K+\delta I_n\)を使うことが普通だ。
そういう訳で数値的不安定性が現れること以外はソフトマージン最大化の場合と2次計画問題を解くコードは全く変わらない。
サポートベクターマシンとPythonによる2次計画問題：導入からスラック変数までに具体的なコードは任せるとしていったんそのカーネル化SVMが定める決定境界の様子を見てみる。
カーネル化SVMが定める決定境界
以上の話はよく知られていることだが、カーネル化SVMの自由度が非常に高いということを図的にみよう。
実験用にこういうデータを作った。これに対していろんなカーネル関数を用いたサポートベクターマシンがどういう決定境界を定めるのかを示していく。以下断らなければ\(C=10\)だ。
まずは次の形の多項式カーネルで\(p\)をいろいろ取った場合を見ていく。
\[ k(x,y)=(1+x^\top y)^p\]
次数が大きくなるにしたがってうねうね度が増していっていることがわかる。4次多項式の場合は見えている範囲では線が2つあるように見える。次に指数関数を使う系
\[ \begin{gather}k(x,y)=\exp\left(-\frac{\|x-y\|^2}2\right)\\k(x,y)=\exp\left(-\frac{|x-y|}2\right)\end{gather}\]
上はガウシアンカーネル、下はラプラスカーネルなどと言われる。
右下で丸3つとはなかななってはくれない。まだまだ試したいカーネルはいろいろあるけどこれについては別で確かめよう。
Python+CVXOPTによるSVMの計算コード例
カーネル化まで至ると前ページのソフトマージンの場合などとは異なりsklearnなどが使えるので2次計画を陽に解くことにこだわる必要は無いのだがここでもCVXOPTモジュールを用いたコード例を示しておく。カーネルトリックによって内積の計算は\(k(x,y)\)の値を直接計算するだけで良いといったが、プログラムを組む際にはもっと別なカーネル関数の与え方が適切だ。
def Kij(X):
return np.exp(-np.sum((X.reshape(n,p,1)- X.T)**2, axis=1))
def Kj(X, x):
return np.exp(-np.sum((X - x)**2, axis=1))このコードはグラム行列\(K\)とベクトル\(\boldsymbol{k}(x),\;k_j(x)=k(x_j,x)\)をデータ行列\(X\)から直接計算している。プログラムの中でカーネル関数にアクセスするのは常にグラム行列の形を介している。ただ、カーネル法の良く知られた問題点として\(n\)が大きい場合にはその記憶容量がバカでかくなることがあるので、その場合には上記の戻り値のメモリを確保すること自体が出来なくなるり、その場合は低ランク近似などによってグラム行列すら介さずに計算を進める必要がある。ここではそういう\(n\)が大きいケースにはあたらない前提でやっている。何がいいたいかというとこの場合でも
def k(x,y):
return np.exp(-np.sum((x-y)**2))のような関数は不要だということだ。紙の上でわかっていればよくてプログラム上は不要になる。これとリスト内包表記などでグラム行列を作ってるとすごく遅い！
cvxoptを使ってカーネルSVMを解くPythonコード例を示す。
def kSVMsolve(X, Y, C=10):
n, p = X.shape
K = Kij(X) + 1e-3*np.eye(n)
P = opt.spdiag([opt.matrix(K), *np.zeros(n+1)])
q = opt.matrix([*np.zeros(n+1), *(C*np.ones(n))])
G = -opt.sparse([
[
opt.matrix(Y.reshape(-1, 1)*np.c_[K, np.ones((n, 1))]),
opt.spmatrix([], [], [], (n, n+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, abstol=abstol, reltol=reltol)
x = np.array(sol['x']).ravel()
alpha, b, xi = x[n:], x[n], x[n+1:] #解ここの部分は先に書き下した2次計画法の問題を素朴に構成してソルバにブチ込んでいるだけだ。スパースな部分はcvxoptに用意されている関数でなるべくスパースに表現している。デフォルトではメッセージと共に効率よく2次計画問題を解いてくれる。主問題と双対問題のコストが一致してれば最適ということだろう。ここに何やってるかの説明があった。
pcost dcost gap pres dres
0: -6.2420e+03 4.8481e+03 2e+04 1e+01 5e+01
1: -9.2347e+02 -2.7643e+02 5e+03 3e+00 1e+01
2: 4.6509e+02 2.5067e+02 1e+03 6e-01 2e+00
3: 6.0075e+02 4.7300e+02 3e+02 1e-01 6e-01
4: 6.2130e+02 5.3830e+02 2e+02 5e-02 2e-01
5: 6.2767e+02 5.6717e+02 1e+02 3e-02 1e-01
6: 6.1883e+02 5.9005e+02 4e+01 6e-03 2e-02
7: 6.0938e+02 6.0168e+02 8e+00 3e-13 5e-10
8: 6.0619e+02 6.0410e+02 2e+00 3e-13 6e-10
9: 6.0521e+02 6.0484e+02 4e-01 3e-13 1e-09
10: 6.0501e+02 6.0500e+02 2e-02 3e-13 4e-09
11: 6.0500e+02 6.0500e+02 3e-04 3e-13 3e-08
Optimal solution found.プログレス表示がいらない場合はどこかで
opt.solvers.options['show_progress'] = Falseを書いておけばいい。最初キーワード引数show_progress=Falseが聞くのかと思ったらそうではないらしい。
実際のところこのページの内容に関しては決定境界の線を引く部分のコードがカーネルSVMを解くところなんかよりも格段に工夫凝らしポイントが高いだがそれはまた別で。
まとめ
サポートベクターマシンとカーネル法を組み合わせると、人が見て分けられそうな分布は大抵分けることが出来るような性質を持っている。ただしある程度適切なカーネル関数の選択やデータを見ての調整が無いと厳しいものは厳しい。
参考文献
福水健次著「カーネル法入門―正定値カーネルによるデータ解析―」、朝倉書店2016. pp. 47, 3.4「サポートベクターマシンの基本事項」