LoginSignup
124

More than 5 years have passed since last update.

pythonでSVM実装

Last updated at Posted at 2016-04-30

はじめに

pythonでSVM(サポートベクターマシン)を実装しました.
教科書として『はじめてのパターン認識』を使いました.

本記事の構成

  • はじめに
  • SVM
    • マージン最大化
    • 主問題
    • ラグランジュ未定乗数法
    • KKT条件
    • 双対問題
  • 解の求め方
  • pythonでの実装
  • 結果
  • おわりに
  • 追記情報

SVM

SVMとは,最大マージンを実現する2クラス線形識別学習法です.
教科書によると,他クラスの識別は $K-1$ 個の一対多SVMで実現することが多いようです.

マージン最大化

次のようなデータ集合を考えます. データの総数を $N$ とすると,$i = 1,...,N$となります.

\begin{align}
& D_L = \bigl\{(t_i, {\boldsymbol{x}}_i)\bigr\} \\
& t_i = \{-1, +1\}
\end{align}

線形識別関数のマージンを $\kappa$ とすると,全ての学習データで次式が成り立ちます.

|\boldsymbol{w}^T{\boldsymbol{x}}_i + b| \geq \kappa

$\kappa$ で正規化したものを $\boldsymbol{w}$, $b$ と置き直し,教師との積をとることで絶対値をはずします.

t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) \geq 1 \tag{1}

図のようにクラス間マージンは,識別超平面に最も近いデータを $\boldsymbol{w}$ 方向に射影した長さの差として求めることができます.

mergin.png

$t_{i} = +1$ のクラスの中で識別超平面に最も近いデータは,$\boldsymbol{w}$ 方向に射影した長さがそのクラス内で最小となり,
$t_{i} = -1$ のクラスの中で識別超平面に最も近いデータは,$\boldsymbol{w}$ 方向に射影した長さがそのクラス内で最大となります.
クラス間マージン $\rho(\boldsymbol{w}, b)$ は次式で表されます.

\begin{align}
\rho(\boldsymbol{w}, b) &= \min_{\boldsymbol{x} \in C_{y=+1}}\frac{{\boldsymbol{w}}^T\boldsymbol{x}}{||\boldsymbol{w}||} - \max_{\boldsymbol{x} \in C_{y=-1}}\frac{{\boldsymbol{w}}^T\boldsymbol{x}}{||\boldsymbol{w}||} \\
\\
&= \frac{1-b}{||\boldsymbol{w}||} - \frac{-1-b}{||\boldsymbol{w}||} \\
\\
&= \frac{2}{||\boldsymbol{w}||} \tag{2}
\end{align}

式(2)より,$||\boldsymbol{w}||$ を最小化すれば,クラス間マージンは最大化されることが分かります.
また,クラス間マージンは識別超平面に最も近いデータだけで決まります.この識別超平面に最も近いデータのことをサポートベクターと呼びます.

主問題

式(1)の制約条件のもと,式(2)に表されるクラス間マージンを最大化することを考えます.これを主問題とします.

評価関数(最小化)

L_p(\boldsymbol{w}) = \frac{1}{2}{\boldsymbol{w}}^{T}\boldsymbol{w}

不等式制約条件

t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1 \geq 0

ラグランジュ未定乗数法

主問題を解くためにラグランジュ未定乗数法を使います.

ラグランジュの未定乗数法(ラグランジュのみていじょうすうほう、英: method of Lagrange multiplier)とは、束縛条件のもとで最適化を行うための数学(解析学)的な方法である。いくつかの変数に対して、いくつかの関数の値を固定するという束縛条件のもとで、別のある1つの関数の極値を求めるという問題を考える。各束縛条件に対して定数(未定乗数、Lagrange multiplier)を用意し、これらを係数とする線形結合を新しい関数(未定乗数も新たな変数とする)として考えることで、束縛問題を普通の極値問題として解くことができる方法である。
wikipediaより引用

具体例

何を言っているかさっぱり分からないので,具体例を挙げて説明します.
$g(x, y) = x^2 + y^2 - 2 = 0$ の条件下で $f(x, y) = y - x$ の最大化を考えます.
これは半径 $\sqrt{2}$ の円周上の $(x, y)$ を使って,直線の切片を最大化するという問題になります.

rag2.png

図から,円と直線が接するとき $f(x, y)$ が制約条件のもと最大値 $2$ を取ることが分かります.
これをラグランジュ未定乗数法を使って解いてみましょう.

\begin{align}
L(x, y, \lambda) &= f(x, y) + \lambda g(x, y) \\
\\
&= y - x + \lambda(x^2 + y^2 - 2) \tag{3}
\end{align}

式(3)のラグランジュ関数を $x, y, \lambda$ で偏微分し,それぞれイコール $0$ とします.

\begin{align}
& \frac{\partial L}{\partial x} = -1 + 2\lambda x = 0 \\
& \frac{\partial L}{\partial y} = 1 + 2\lambda y = 0 \\
& \frac{\partial L}{\partial \lambda} = x^2 + y^2 - 2 = 0
\end{align}

1,2番目の式を3番目の式に代入します.

\begin{align}
& \Bigl(\frac{1}{2\lambda}\Bigr)^2 + \Bigl(-\frac{1}{2\lambda}\Bigr)^2 - 2 = 0 \\
\\
& \lambda = \pm \frac{1}{2}, x = \pm 1, y = \mp 1
\end{align}

よって,$g(x, y) = 0$ の制約条件下で $f(x, y) = y - x$ は最大値 $2$ を取ります.
この最大値について3次元で考えてみます.3次元空間では,$z = f(x, y) = y - x$ は斜面に,$g(x, y) = x^2 + y^2 - 1 = 0$ は円柱になります.
$g(x, y) = 0$ の制約条件下で $z = y - x$ の最大値を求めたいので,斜面と円柱が交わる面に関して,$z$ 座標の最大値を見つければ良いこととなります.
つまり,$x^2 + y^2 - 2 = 0$ の円を $z = y - x$ 平面に射影したときの $z$ 軸方向の最大値が求めたい値です.
これを図で表します.真横から見たとき,斜面と円柱が交わっている部分の中で,$z$ 軸方向の最大値は $2$ であることが分かり,これは計算結果と一致します.

lagrange.png

解説

なぜラグランジュ未定乗数法で,最小値を求められるのか解説してみます.(誤りを含むかもしれません)
制約条件 $g(\boldsymbol{x}) = 0$(以下,制約面とする)のもと,評価関数 $f(\boldsymbol{x})$ が $\boldsymbol{x}^\star$ において最大値をとるとき,$\boldsymbol{\nabla}f(\boldsymbol{x})$ は制約面に対して垂直でなければなりません.
なぜなら,垂直でないと仮定すると,制約面の方向に沿って $f(\boldsymbol{x})$ をさらに大きくすることができるからです.また,$\boldsymbol{\nabla}g(\boldsymbol{x})$ も制約面に対して垂直になります.
よって $\boldsymbol{x}^\star$ において $\boldsymbol{\nabla}f, \boldsymbol{\nabla}g$ は平行なベクトルとなり,次式が成り立ちます.

\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0 \tag{4}

ここで,ラグランジュ関数を以下のように定義します.

L(\boldsymbol{x}, \lambda) \equiv f(\boldsymbol{x}) + \lambda g(\boldsymbol{x})

ラグランジュ関数 $L(\boldsymbol{x}, \lambda)$ を $\boldsymbol{x}, \lambda$ で偏微分することで式(4)および制約条件を導くことができます.
これを解くことで停留点 $\boldsymbol{x}^\star$ が得られます.この $\boldsymbol{x}^\star$ の中に,$f(\boldsymbol{x})$ を最大化する $\boldsymbol{x}$ が存在します.

KKT条件

では,このラグランジュの未定乗数法を主問題に適用します.
ラグランジュ関数を次式のように定義します.

\tilde{L}_{p}(\boldsymbol{w}, b, \boldsymbol{\alpha}) = \frac{1}{2}{\boldsymbol{w}}^T\boldsymbol{w} - \sum_{i=1}^{N}\alpha_i(t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1) \tag{5}

ここで,$\boldsymbol{\alpha} = {(\alpha_1, ..., \alpha_N)}^T (\alpha_i \geq 0)$ であり,$\alpha_i$ がラグランジュ未定乗数です.
この最適化問題の解 $\boldsymbol{w}_0, b_0$ は,以下のKKT条件を満たす解として得られます.

< KKT条件 >

\begin{align}
& \left.\frac{\partial \tilde{L}_p(\boldsymbol{w}, b, \boldsymbol{\alpha})}{\partial \boldsymbol{w}}\right|_{\boldsymbol{w} = \boldsymbol{w}_0} = \boldsymbol{w}_0 - \sum_{i=1}^{N}\alpha_{i}t_i{\boldsymbol{x}}_i = 0 \tag{6} \\
& \frac{\partial \tilde{L}_p(\boldsymbol{w}, b, \boldsymbol{\alpha})}{\partial b} = \sum_{i=1}^{N}\alpha_{i}t_i = 0 \tag{7} \\
& \\
& t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1 \geq 0 \tag{8} \\
& \\
& \alpha_i \geq 0 \tag{9} \\
& \\
& \alpha_i(t_i(\boldsymbol{w}^T{\boldsymbol{x}}_i + b) - 1) = 0 \tag{10}
\end{align}

式(6),(7)はラグランジュ関数の偏微分イコール $0$,式(8)は主問題の制約条件,式(10)は相補性条件です.
式(10)の相補性条件より

t_i(\boldsymbol{w}^T\boldsymbol{x}_i + b) - 1 > 0 \Rightarrow \alpha_i = 0 \\
t_i(\boldsymbol{w}^T\boldsymbol{x}_i + b) - 1 = 0 \Rightarrow \alpha_i \geq 0

相補性条件は,サポートベクター以外のデータによる制約条件の効果はなく、サポートベクターだけで識別超平面が決定されることを示しています.
(調べてみるとKKT条件はかなり奥が深そうなので,これくらいのふんわりとした理解で次に進むことにします)

双対問題

式(5)を展開し,最適解および式(6),(7)を代入します.

\begin{align}
L_{d}(\boldsymbol{\alpha}) &= \frac{1}{2}{\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} - \sum_{i=1}^{N}\alpha_{i}t_{i}{\boldsymbol{w}_{0}}^{T}{\boldsymbol{x}}_{i} - b\sum_{i=1}^{N}\alpha_{i}t_{i} + \sum_{i=1}^{N}\alpha_{i} \\
&= \frac{1}{2}{\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} - {\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} - 0 + \sum_{i=1}^{N}\alpha_{i} \\
&= \sum_{i=1}^{N}\alpha_{i} - \frac{1}{2}{\boldsymbol{w}_{0}}^{T}\boldsymbol{w}_{0} \\
&= \sum_{i=1}^{N}\alpha_{i} - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}t_{i}t_{j}{\boldsymbol{x}_i}^T{\boldsymbol{x}}_j
\end{align}

式(6)より,最適解がラグランジュ未定乗数と学習データの線形和で表されることから,$\alpha_{i}$ を求める問題に置き換えることができます.
最適なラグランジュ未定乗数は,$L_{d}(\boldsymbol{\alpha})$ を最大にする $\boldsymbol{\alpha}$ により得られます.これを主問題に対する双対問題と呼びます.
$N$ 個の $1$ を並べたベクトルを $\boldsymbol{1}$,学習データで作られた行列を $\boldsymbol{H}$,教師データベクトルを $\boldsymbol{t}$ とします.

\begin{align}
& \boldsymbol{1} = (1,...,1)^T \\
\\
& \boldsymbol{H} = (H_{ij} = {t_it_j\boldsymbol{x}_i}^T\boldsymbol{x}_j) \\
\\
& \boldsymbol{t} = (t_1,...,t_N)^T
\end{align}

評価関数(最大化)

L_{d}(\boldsymbol{\alpha}) = \boldsymbol{\alpha}^T\boldsymbol{1} - \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{H}\boldsymbol{\alpha}

等式制約条件

\boldsymbol{\alpha}^T\boldsymbol{t} = 0

双対問題を解いて得られたラングランジュ未定乗数の最適解を $\tilde{\boldsymbol{\alpha}} = (\tilde{\alpha}_1,...,\tilde{\alpha}_N)^T$ とすると,最適解は次のように求められます.

\begin{align}
& \boldsymbol{w}_0 = \sum_{i=1}^{N}\tilde{\alpha}_{i}t_i{\boldsymbol{x}}_i \\
& b_0 = \frac{1}{t_s} - {\boldsymbol{w}_0}^T\boldsymbol{x}_s
\end{align}

$\boldsymbol{x}_s, t_s$ はサポートベクターのデータおよび教師を表します.

解の求め方

主問題を双対問題に置き換え,最適な $\tilde{\boldsymbol{\alpha}}$ を求めれば,マージン最大化は実現できることを導きました.
$L_{d}(\boldsymbol{\alpha})$ を最大にする $\boldsymbol{\alpha}$ をプログラムで求めるために最急降下法を使います.
双対問題の制約条件を満たすよう罰金項を加え,以下のように $L_{d}(\boldsymbol{\alpha})$ を再定義しました.

L_{d}(\boldsymbol{\alpha}) = \boldsymbol{\alpha}^T\boldsymbol{1} - \frac{1}{2}\boldsymbol{\alpha}^T\boldsymbol{H}\boldsymbol{\alpha} - \frac{1}{2}\beta\boldsymbol{\alpha}^T\boldsymbol{t}\boldsymbol{t}^T\boldsymbol{\alpha}

では,実際に最急降下法で $\alpha_i$ を更新します.

\begin{align}
{\alpha_i}' &= {\alpha_i} + \eta\frac{L_{d}(\boldsymbol{\alpha})}{\alpha_i} \\
&= {\alpha_i} + \eta(1 - \sum_{j=1}^{N}\alpha_{j}t_{i}t_{j}{\boldsymbol{x}_i}^T{\boldsymbol{x}}_j - \beta\sum_{j=1}^{N}\alpha_{j}t_{i}t_j)
\end{align}

ここまで難しい計算をしてきましたが,最終的に得られる数式は美しいですね.
これで $L_{d}(\boldsymbol{\alpha})$ を最大にする $\boldsymbol{\alpha}$ は得られます.

pythonでの実装

以下のように実装しました.
$\alpha_i$ の更新が簡単にできます.また,罰金項の効果で,$\boldsymbol{\alpha}^T\boldsymbol{t} = 0$ の条件を満たすように更新が行われます.

svm.py
import numpy
from matplotlib import pyplot
import sys

def f(x, y):
    return x - y

if __name__ == '__main__':

    param = sys.argv

    numpy.random.seed()
    N = 30
    d = 2
    X = numpy.random.randn(N, d)
    T = numpy.array([1 if f(x, y) > 0 else - 1 for x, y in X])
    alpha = numpy.zeros(N)
    beta = 1.0
    eta_al = 0.0001 # update ratio of alpha
    eta_be = 0.1 # update ratio of beta
    itr = 1000

    for _itr in range(itr):
        for i in range(N):
            delta = 1 - (T[i] * X[i]).dot(alpha * T * X.T).sum() - beta * T[i] * alpha.dot(T)
            alpha[i] += eta_al * delta
        for i in range(N):
            beta += eta_be * alpha.dot(T) ** 2 / 2

    index = alpha > 0
    w = (alpha * T).T.dot(X)
    b = (T[index] - X[index].dot(w)).mean()

    if '-d' in param or '-s' in param:
        seq = numpy.arange(-3, 3, 0.02)
        pyplot.figure(figsize = (6, 6))
        pyplot.xlim(-3, 3)
        pyplot.ylim(-3, 3)
        pyplot.plot(seq, -(w[0] * seq + b) / w[1], 'k-')
        pyplot.plot(X[T ==  1,0], X[T ==  1,1], 'ro')
        pyplot.plot(X[T == -1,0], X[T == -1,1], 'bo')

        if '-s' in param:
            pyplot.savefig('graph.png')

        if '-d' in param:
            pyplot.show()

結果

下の図のような結果が得られました.データ点の色がクラスを表しています.
計算により求められた直線を一緒に描画したところ,それっぽい結果は得られました.

result.png

おわりに

不等式条件下でマージン最大化を目的とした主問題を設定し,これをラグランジュ未定乗数法を使うことで双対問題へ置き換え,最適な $\boldsymbol{w}_0, b_0$ を求められることが分かりました.

追記情報

  • 2016/05/01 @skyshk さんからコメントをいただいてコードを修正しました.

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
124