サポートベクターマシンの実装（その１）

ラグランジュ係数を2次計画問題で解くところが肝。
SMOアルゴリズムとかあるけど、
とりあえず最急降下法を試す。

やることは、ラグランジュ関数
$L(\alpha) = \sum_{i=1}^n \alpha_i - \frac{1}{2} \sum \alpha_i \alpha_j y_i y_j \mathbf{{x_i}^t} \mathbf{x_j}$
の、制約条件
$\sum_{i=1}^n \alpha_i y_i = 0$
$\alpha_i \geq 0$
のもとでの最大化である。

更新式は、
$\alpha_i ' = \alpha_i - \rho \frac{\partial L}{\partial \alpha_i}$ 　より
$\alpha_i ' = \alpha_i - \rho (1- \sum_{j=1}^n \alpha_j y_i y_j \mathbf{{x_i}^t} \mathbf{x_j} )$

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
from pylab import *

#　サンプルデータの生成（２次元正規分布）
mean1 = [0,0]; cov1 = [[4,0],[0,100]]; N1 = 100
X1 = np.random.multivariate_normal(mean1,cov1,N1)
mean2 = [30,-30]; cov2 = [[1,20],[20,50]]; N2 = 100
X2 = np.random.multivariate_normal(mean2,cov2,N2)
X = np.concatenate((X1,X2))
# 描画
x,y = X.T
plt.plot(x,y,'k.'); plt.axis('equal'); #plt.show()

#ラグランジュ係数の解ベクトル
alpha = []
for k in range(100):
    alpha.append(0.01)
#print alpha
A = np.concatenate((alpha,alpha))

y1 = []
for k in range(100):
    y1.append(1)
y2 = []
for k in range(100):
    y2.append(-1)
Y = np.concatenate((y1,y2))

#最急降下法(勾配上昇法)でラグランジュ関数を最大化する
rho = 0.000001

for n in range(10):
    for i in range(200):
        sum = 0
        for j in range(200):
            sum += Y[i]*A[j]*Y[j]*(np.dot(X[i],X[j]))
        A[i] = A[i] - rho * (1 - sum)
        print i
print A

計算量が多いので特徴ベクトルの個数＝各クラス１００個、
繰り返し回数10回にしてみた。

サポートベクトルらしきものが見つからないので保留。