Python
機械学習
svm

pythonでSVM実装

More than 1 year has passed since last update.

はじめに

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 さんからコメントをいただいてコードを修正しました.