SVMを実装してみた

授業でSVMについて習ったけど、実際に実装したことなかったからやってみた。簡単って言われてるけど、制約付き2次計画問題の実装が結構大変だった（収束しないケースとかたくさんあったり、制約条件を遵守したり）

参考にしたのは、以下の本やらページやら
http://www.amazon.co.jp/%E3%82%B5%E3%83%9D%E3%83%BC%E3%83%88%E3%83%99%E3%82%AF%E3%82%BF%E3%83%BC%E3%83%9E%E3%82%B7%E3%83%B3%E5%85%A5%E9%96%80-%E3%83%8D%E3%83%AD-%E3%82%AF%E3%83%AA%E3%82%B9%E3%83%86%E3%82%A3%E3%82%A2%E3%83%8B%E3%83%BC%E3%83%8B/dp/4320121341/ref=sr_1_1?ie=UTF8&qid=1354439820&sr=8-1
http://home.hiroshima-u.ac.jp/tkurita/lecture/svm/index.html

制約条件の中でも厄介なのが、ラベルのベクトルyと重み係数ベクトルαの内積が0になる必要があるという制約で、普通にやると結構この制約を満たすだけでも一苦労…
そこで考えたのが、最初からyに直交する空間上で最適化するという案。具体的に基底を定めてやると、意外と簡単にαからβという別の重み係数の最適化問題に落とせたので、それで実装した。ちなみに最適化エンジンは単なる勾配法。
完成はしたけど学習率とかカーネル等を上手く調整しないと収束しない場合が結構ある。理由はよく分かってないけど、経験的にはサポートベクターの数が多いと収束しやすい傾向があると予想してる（性質としては良くないけどw）

コードはこんな感じになった。サポートベクター表示させてなかったり、サポートベクターじゃない奴も含めてクラス分類させてたりあるけど…
行列による一斉更新ってところは、βを逐次的に変えるより速いんだけど収束しなかったり収束が遅かったりするからコメントアウトしてある。

# -*- coding: utf-8 -*-

from pylab import *
from numpy.linalg import *

LABEL = 0
VECTOR = 1

# y = func(x)が境界線
def make_input_data(Nsample, func):
    data = []

    #正例のデータ作成
    while len(data) < Nsample / 2:
        xs = rand(2)
        if xs[1] >= func(xs[0]):
            data.append((1, xs))
            
    #負例のデータ作成
    while len(data) < Nsample:
        xs = rand(2)
        if xs[1] < func(xs[0]):
            data.append((-1, xs))

    return data

#係数alphaの最適化を行う
#要件: 最後のサンプルのラベルが-1であること
def optimize(samples, kernel):

    #bの次元数
    dim = len(samples) - 1

    #alphaの代わりにbを最適化する
    #最適化問題: min f(b) = b^TGb - c^Tb
    
    c = mat(zeros((dim, 1)))
    for i in xrange(dim):
        c[i] = 1 + samples[i][LABEL]

    #kernel(x_i, x_n)を求める
    gn = zeros(dim + 1)
    for i in xrange(dim + 1):
        gn[i] = kernel(samples[i][VECTOR], samples[dim][VECTOR])
    
    G = mat(zeros((dim, dim)))
    for y in xrange(dim):
        for x in xrange(y, dim):
            G[y, x] = G[x, y] = kernel(samples[x][VECTOR], samples[y][VECTOR]) \
                - gn[x] - gn[y] + gn[dim]

    #bの初期化
    b = mat(zeros((dim, 1)))
    for i in xrange(dim):
        b[i] = 0.5 * samples[i][LABEL]
    
    before_b = mat(ones((dim, 1)))
    dist_eps = 1.0e-4
    nu = 1.0e-4
    iteration = 0
    max_iteration = 1e5

    while norm(b - before_b) > dist_eps and iteration < max_iteration:

        if iteration % 1000 == 0:
            print norm(b - before_b), iteration

        before_b = b.copy()

        #bを行列による一斉更新
        #d = dot(G, b) - c
        #b -= nu * d

        #bを逐次更新
        for i in range(dim):
            d = dot(G[i,:], b) - c[i]
            b[i] -= nu * d
        
        #制約による修正 (bi * yi >= 0)
        for i in range(dim):
            if samples[i][LABEL] * b[i] < 0:
                b[i] = 0

        #制約による修正 (sum(b) >= 0)
        sum_b = sum(b)
        if sum_b < 0:
            Np = len(filter(lambda item: item[LABEL] == 1, samples))
            for i in xrange(dim):
                if samples[i][LABEL] == 1:
                    b[i] -= sum_b / Np
                    
        iteration += 1

    print iteration
        
    #biをalpha_iに直す
    alphas = zeros(dim + 1)
    for i in xrange(dim):
        alphas[i] = samples[i][LABEL] * b[i]
    alphas[dim] = sum(b)

    #alpha_iの最大値より十分小さい要素は0
    alpha_threashold = max(alphas) * 1.0e-5
    for i in xrange(len(alphas)):
        if alphas[i] < alpha_threashold:
            alphas[i] = 0.0
            
    return alphas

#最適化したalphaを用いてサンプルの最適化を行う
def make_classifier(kernel, alphas, samples):

    def weighted_sum_kernel(x):
        ans = 0.0
        for i in range(len(samples)):
            ans += kernel(samples[i][VECTOR], x) * alphas[i] * samples[i][LABEL]
        return ans
        
    psamples = filter(lambda sample: sample[LABEL] == 1, samples)
    min_psection = min(map(lambda sample: weighted_sum_kernel(sample[VECTOR]), psamples))
    nsamples = filter(lambda sample: sample[LABEL] == -1, samples)
    max_nsection = max(map(lambda sample: weighted_sum_kernel(sample[VECTOR]), nsamples))
    section = - (max_nsection + min_psection) / 2.0
    
    def classify(x):
        return weighted_sum_kernel(x) + section

    return classify    

#カーネルの一例
def gauss_kernel(sigma):
    def kernel(x1, x2):
        return exp(-(norm(x1 - x2)) / ((2 * sigma) ** 2))
    
    return kernel

def inner_product():
    def kernel(x1, x2):
        return dot(x1.T, x2)
    return kernel

#境界面の関数の一例
def target(x):
    return sin(4 * x) / 2.0 + 0.25

def target2(x):
    return 1 - x

#実際の使用例
Nsample = 200

func = target

samples = make_input_data(Nsample, func)
kernel = gauss_kernel(0.5)
alphas = optimize(samples, kernel)

classify = make_classifier(kernel, alphas, samples)

for i in range(len(samples)):
    print samples[i][LABEL], classify(samples[i][VECTOR]), alphas[i]

class1 = [item[VECTOR] for item in samples if item[LABEL] == 1]
class2 = [item[VECTOR] for item in samples if item[LABEL] == -1]

xs1 = [v[0] for v in class1]
ys1 = [v[1] for v in class1]
xs2 = [v[0] for v in class2]
ys2 = [v[1] for v in class2]

plot(xs1, ys1, 'ro')
plot(xs2, ys2, 'bo')

xs = arange(0, 1, 0.01)

size = len(xs)

mesh_z = zeros((size, size))
xy = zeros(2)

for y in range(size):
    for x in range(size):
        xy[0], xy[1] = xs[x], xs[y]
        mesh_z[y, x] = classify(xy)

mesh_x, mesh_y = meshgrid(xs, xs)
contourf(mesh_x, mesh_y, mesh_z, 0)
show()

これが結果の一例で、上手くいけばこうなった。

SVMはやっぱすごいけど、実際使ったり実験するならやっぱパッケージ使うよなーと実感した