サポートベクターマシン（SVM）について書きたいと思います。最初にサポートベクターマシンのざっくりとした理論について述べて、そのあとに実装について書きたいと思います。

線形SVM

与えられた線形分離可能なデータを2つのクラスに分けることを考えます。このとき、識別平面の引き方は割と遊びがあります。例えば、以下のようなデータを2クラスに分けるときには

f:id:kanoeuma_310mo:20190922142757p:plain:w300 — こんな引き方から

f:id:kanoeuma_310mo:20190922142819p:plain:w300 — こんな引き方まで

色々な引き方が考えられます。この中で、最適な引き方を考えます。

マージンという距離を定義します。これは2つのクラスのうち、異なるクラスに最も近いデータ同士の距離のことみたいなイメージです。このマージンを最大にするような識別平面を定義することで最適な引き方とします。

f:id:kanoeuma_310mo:20190922144704p:plain:w300 — こんな感じ

マージンの定義とか識別平面とかについて

重みwとバイアスbを用いて識別平面を以下のように定義します。

$w \cdot x - b = 0$

このとき、マージンを以下の2つの平面の距離とします。この2つの平面は、分類されるクラスの端っこの点を通る平面を表しています。

$w \cdot x - b = 1 , w \cdot x - b = -1$

すなわち、マージンは $\frac{2}{||w||}$ で表されます。（点と平面の距離とか久しく使ってないので忘れてた）

よって以下の条件を満たしながらマージンを最大化（要するに $||w||$ を最小化）すれば良いことになります。yはデータのクラスラベルです。（はてなブログでのcasesの使い方がいまいちわからん）

$w \cdot x - b \geq 1(y=1), w \cdot x - b \leq -1(y=-1)$

これは $y(w \cdot x - b) \geq 1$ にまとめられるので、結果として $y(w \cdot x - b) \geq 1$ という条件下で $||w||$ を最小化する問題とすることができます。

この問題はラグランジュ未定乗数 $\alpha_i$ によって以下のような双対形式に変形して解くことができます。（ここでは $||w||$ ではなく $\frac{1}{2}||w||^2$ を最小化することを考えます）ここの変形の妥当性をあんまりしっかり理解していないのごめんなさい...

$L = \frac{1}{2}||w||^2 - \sum_{i=1}^n\alpha_i(y_i(w \cdot x_i -b)-1)$ （このLを $\alpha$ について最大化し、 $w, b$ について最小化する）

また、この式には以下の条件が付随します。

$\frac{\partial L}{\partial w}=0, \frac{\partial L}{\partial b}=0, \alpha_i \geq 0, \alpha_i(y_i(w_i \cdot x_i - b) - 1) = 0$

$\frac{\partial L}{\partial w}=0$ より $w = \sum \alpha_i y_i x_i$ , $\frac{\partial L}{\partial b}=0$ より $\sum \alpha_i y_i=0$ が成立する必要があります。

したがって、元々の問題は以下のLを最大化する変形することができます。（はてなブログってalign環境使えないんです？）

$L = \frac{1}{2}||w||^2 - \sum_{i=1}^n w \cdot \alpha_i y_i x_i$

$- \sum_{i=1}^n \alpha_i y_i b$

$+ \sum_{I=1}^n \alpha_i$

$= \sum_{I=1}^n \alpha_i - \frac{1}{2} \sum \alpha_i \alpha_j y_i y_j x_i^{T} x_j$

subject to $\alpha_i \geq 0, \sum_{I=1}^n \alpha_i y_i = 0$

これを二次計画問題として解きます。（二次計画問題を解くところはライブラリにぶん投げたので省きます。理解してないし）

プログラムによる実装

要するに二次計画問題にぶっ込めばいいということになります。cvxoptを使っています。

線形SVMを二次計画問題として解釈してcvxoptに突っ込むとこんな感じになります。

h=x*l
qpP = cvxopt.matrix(h.T.dot(h))
qpq = cvxopt.matrix(-np.ones(n), (n, 1))
qpG = cvxopt.matrix(-np.eye(n))
qph = cvxopt.matrix(np.zeros(n), (n, 1))
qpA = cvxopt.matrix(l.astype(float), (1, n))
qpb = cvxopt.matrix(0.)
cvxopt.solvers.options['abstol'] = 1e-5
cvxopt.solvers.options['reltol'] = 1e-10
cvxopt.solvers.options['show_progress'] = False
res=cvxopt.solvers.qp(qpP, qpq, qpG, qph, qpA, qpb)
alpha = np.reshape(np.array(res['x']),-1)

xはデータ、lはラベル、nはデータ数です。

ソフトマージン

線形SVMは与えられたデータが線形分離可能という割ときつい制約がありました。これをある程度解決するのがソフトマージンです。

ソフトマージンは $y(w \cdot x -b) \geq 1$ を満たせないデータを許容する代わりに各データにペナルティとしてスラック変数 $\zeta_i$ を以下のように与えます。

$y_i(w \cdot x_i -b) \geq 1 - \zeta_i$

$\zeta_i$ が0のときは線形SVMと同様にマージン内で正しく識別できていることを表します。0以上1以下の場合はマージンの境界を超えてしまうものの、識別平面を踏み越えているわけではないので（識別平面は=0の式）正しく識別できていることを表します。1以上になってしまった場合は識別平面を踏み越えて誤識別されていることを表します。

線形SVMのときには $\frac{1}{2}||w||^2$ を最小にすることを考えましたが、スラック変数を導入したことで最小化する関数は以下のように変化します。

$\frac{1}{2}||w||^2 + C \sum\zeta_i$

subject to $y_i(w \cdot x_i -b) \geq 1 - \zeta_i, \zeta_i \geq 0$

Cはペナルティである $\zeta_i$ をどれだけ重く見るかを表す定数です。重く見れば見るほど線形SVMライクになるということであってると思います。

これを線形SVMのときと同じようにラグランジュ未定乗数とかでごにょごにょすると、以下の条件でLを最大化する問題に変形することができます。

$L = \sum \alpha_i - \frac{1}{2} \sum \alpha_i \alpha_j y_i y_j x_i^{T} x_j$

subject to $0 \leq \alpha_i \leq C, \sum_{I=1}^n \alpha_i y_i = 0$

実際のところ、ソフトマージンを導入したことによる差は $\alpha_i \leq C$ の条件が追加されることになります。

プログラムによる実装

線形SVMのときの二次計画問題の条件に $\alpha_i \leq C$ が追加されるだけなので以下のように計算すればいいです。

h=x*l
C = 1.0

qpP = cvxopt.matrix(h.T.dot(h))
qpq = cvxopt.matrix(-np.ones(n), (n, 1))
qpG = cvxopt.matrix(np.vstack((-np.eye(n), np.eye(n))))
qph = cvxopt.matrix(np.hstack((np.zeros(n), C*np.ones(n))), (2*n, 1))
qpA = cvxopt.matrix(l.astype(float), (1, n))
qpb = cvxopt.matrix(0.)
cvxopt.solvers.options['abstol'] = 1e-5
cvxopt.solvers.options['reltol'] = 1e-10
cvxopt.solvers.options['show_progress'] = False
res=cvxopt.solvers.qp(qpP, qpq, qpG, qph, qpA, qpb)
alpha = np.reshape(np.array(res['x']),-1)

xはデータ、lはラベル、nはデータ数です。不等式の条件の部分で配列を結合しています。

カーネルトリック

ソフトマージンとかではもうどうしようもないほど線形分離不可能だけど分離自体は可能である場合に対応する方法を考えます。

その場合に使うのがガウスカーネル $k(x, y) = \exp(-||x-y||^2/2$ sigma²)などのカーネル関数です。（はてなブログのtex記法なんか文字数多くなると機能しなくなる気がする？）これを今まで内積で計算していた部分に代わりに代入することで線形分離不可能なデータにも対応できます。

データをxを単純に非線型変換してから内積を取るとカーネル関数が表す関数への非線形変換は一般的には難しいです。でも、内積をカーネルごと入れ替えることでxの非線形変換に落とし込むことなく変換をすることができます。

プログラムによる実装

今まで内積を適用していた部分をガウスカーネルとかに変えてやれば良いです。カーネル関数はガウスカーネルを使っています。

h=x*l
C = 1.0
sigma = 1.0



def gauss_kernel(x, y):
    return np.exp(-np.linalg.norm(x-y)**2 / sigma)

P = np.zeros((n, n))

for i in range(n):
    for j in range(n):
        P[i][j] = l[i]*l[j]*gauss_kernel(x.T[i], x.T[j])

qpP = cvxopt.matrix(P)
qpq = cvxopt.matrix(-np.ones(n), (n, 1))
qpG = cvxopt.matrix(np.vstack((-np.eye(n), np.eye(n))))
qph = cvxopt.matrix(np.hstack((np.zeros(n), C*np.ones(n))), (2*n, 1))
qpA = cvxopt.matrix(l.astype(float), (1, n))
qpb = cvxopt.matrix(0.)
cvxopt.solvers.options['abstol'] = 1e-5
cvxopt.solvers.options['reltol'] = 1e-10
cvxopt.solvers.options['show_progress'] = False
res=cvxopt.solvers.qp(qpP, qpq, qpG, qph, qpA, qpb)
alpha = np.reshape(np.array(res['x']),-1)

まとめ

以上です。説明省いたり、実装とかほざいて二次計画問題のライブラリを使ってるだけだったりとグダグダですがある程度まとめられたかなとは思います。

庚午里藻の日記

見た映画とかアニメの備忘録にしたり、パソコンいじったことのメモにしたり

サポートベクターマシンについて

線形SVM

マージンの定義とか識別平面とかについて

プログラムによる実装

ソフトマージン

プログラムによる実装

カーネルトリック

プログラムによる実装

まとめ