ソフトマージンSVMを書いてみた

numpyの練習を兼ねてソフトマージンSVMを書きました．
行列演算に関してはRと比べても大きな不自由は感じないですね．

最適化部分は二次計画法のソルバー(openopt, cvxopt)に丸投げ．
openoptから呼べるソルバーのうちデフォルトで含まれているqlcpではデータによっては解けなかったのでcvxoptを利用しました．

ソフトマージンSVM - 人工知能に関する断創録の記事が非常にわかりやすくておすすめです．

ソースコード

適当にガウス分布を2つ作って学習してます．
必要なライブラリ，numpy，openopt，cvxopt．

import numpy as np
import openopt as oo

#settings
C = 1
EPS = 0.00001

#data generation
x1 = np.random.multivariate_normal([0, 0], np.eye(2), 10)
x2 = np.random.multivariate_normal([2, 2], np.diag([2, 2]), 10)
x = np.vstack((x1, x2))
t = np.array([-1] * x1.shape[0] + [1] * x2.shape[0])

N = x.shape[0]

#optimize
qp = oo.QP((t * np.dot(x, x.T)).T * t, -np.ones(N), Aeq = t, beq = 0, lb = np.zeros(N), ub = np.array([C] * N))
result = qp._solve('cvxopt_qp', iprint = 0)
a = result.xf

#parameter estimation
w = (x.T * t * a).sum(axis = 1)
support_vector = a > EPS
b = (sum(t[support_vector]) - sum(np.dot(w, x[support_vector].T))) / sum(support_vector)

#output
print np.dot(w, x.T) + b