この記事は何?
サポートベクターマシン（Support Vector Machine: SVM）では以下の目的関数を最小化します．
\min_{\bf w}\frac{\lambda}{2}\|{\bf w}\|^2+\frac{1}{m}\sum_{({\bf x},y)\in S}l({\bf w}; ({\bf x}, y)) \\
l({\bf w}; ({\bf x},y))=\max(0,1-y\left<{\bf w},{\bf x} \right>)
オンライン機械学習に，この目的関数を確率的勾配降下法で最小化するアルゴリズムが載っていましたが，カーネルトリックには対応していませんでした．
そこで，確率的勾配降下法によるカーネルトリックを用いたSVMの実装を探していたところ，Pegasos: Primal Estimated sub-GrAdient SOlver for SVMという論文を見つけました．
この記事では，論文中で示されていたアルゴリズムを実装してみます．
アルゴリズムの入力と出力
論文中のFig.3を実装します．
アルゴリズムの入力は以下の通りです．
- 訓練データ: $S=\lbrace(x_i,y_i)\rbrace_{i=1}^m$
- 学習率: $\lambda$（正確には，$\eta=\frac{1}{\lambda t}$という学習率が使われます．）
- エポック数: $T$
- カーネル関数: $K(x_i,x_j)$
Fig.3中のChoose $i_t \in \lbrace 0,\dots,|S| \rbrace$の部分はChoose$i_t \in \lbrace 1,\dots,|S|\rbrace$の間違いだと思います（たぶん）．
アルゴリズムの出力は，$\alpha_{T+1}$という要素数$m$のベクトルです．
得られた$\alpha_{T+1}$を用いることで，以下のように予測を行うことができます．
\hat{y}={\bf w}_{T+1}^{\rm T}\phi({\bf x})=\frac{1}{\lambda t}\sum_{j=1}^{m}\alpha_{T+1}[j]y_jK({\bf x}_j,{\bf x})
pythonによる実装例
sklearnライクに実装してみました．また，バイアス項を加えるようにアルゴリズムを修正しています．
ソースコードはここにアップロードしてあります．
import numpy as np
class SGDSVC(object):
def __init__(self,
kernel="rbf", lmd=1e-1, gamma=0.1, bias=1.0, max_iter=100):
if kernel not in self.__kernel_dict:
print(kernel + " kernel does not exist!\nUse rbf kernel.")
kernel = "rbf"
if kernel == "rbf":
def kernel_func(x, y):
return self.__kernel_dict[kernel](x, y, gamma=gamma)
else:
kernel_func = self.__kernel_dict[kernel]
self.kernel = kernel_func
self.lmd = lmd
self.bias = bias
self.max_iter = max_iter
def __linear_kernel(x, y):
return np.dot(x, y)
def __gaussian_kernel(x, y, gamma):
diff = x - y
return np.exp(-gamma * np.dot(diff, diff))
__kernel_dict = {"linear": __linear_kernel, "rbf": __gaussian_kernel}
def fit(self, X, y):
def update_alpha(alpha, t):
data_size, feature_size = np.shape(self.X_with_bias)
new_alpha = np.copy(alpha)
it = np.random.randint(low=0, high=data_size)
x_it = self.X_with_bias[it]
y_it = self.y[it]
if (y_it *
(1. / (self.lmd * t)) *
sum([alpha_j * y_it * self.kernel(x_it, x_j)
for x_j, alpha_j in zip(self.X_with_bias, alpha)])) < 1.:
new_alpha[it] += 1
return new_alpha
self.X_with_bias = np.c_[X, np.ones((np.shape(X)[0])) * self.bias]
self.y = y
alpha = np.zeros((np.shape(self.X_with_bias)[0], 1))
for t in range(1, self.max_iter + 1):
alpha = update_alpha(alpha, t)
self.alpha = alpha
def decision_function(self, X):
X_with_bias = np.c_[X, np.ones((np.shape(X)[0])) * self.bias]
y_score = [(1. / (self.lmd * self.max_iter)) *
sum([alpha_j * y_j * self.kernel(x_j, x)
for (x_j, y_j, alpha_j) in zip(
self.X_with_bias, self.y, self.alpha)])
for x in X_with_bias]
return np.array(y_score)
def predict(self, X):
y_score = self.decision_function(X)
y_predict = map(lambda s: 1 if s >= 0. else -1, y_score)
return y_predict
実行例
irisデータセットに対する予測例を示します．
実装の正しさを確認するために，訓練時に用いたデータをそのままテストに用います．
パラメータは$\lambda=0.1, T=100$で，カーネル関数にはガウスカーネル（$\gamma=0.1$）を使いました．
また，訓練データは平均が0，分散が1になるようなスケーリングを行っています．
# -*- coding: utf-8 -*-
import SGDSVC
from sklearn import datasets
from sklearn import preprocessing
from sklearn import metrics
if __name__ == '__main__':
iris = datasets.load_iris()
X = iris.data[:100]
y = iris.target[:100]
# ラベルは1, -1
y = map(lambda l: 1 if l == 1 else -1, y)
# スケーリング
scaler = preprocessing.StandardScaler()
X = scaler.fit_transform(X)
clf = SGDSVC.SGDSVC()
clf.fit(X, y)
y_pred = clf.predict(X)
print("Confusion Matrix")
print(metrics.confusion_matrix(y, y_pred))
print(metrics.classification_report(y, y_pred))
実行結果は以下の通りです．
Confusion Matrix
[[50 0]
[ 0 50]]
precision recall f1-score support
-1 1.00 1.00 1.00 50
1 1.00 1.00 1.00 50
avg / total 1.00 1.00 1.00 100
上で述べたように，訓練時に用いたデータをそのままテストに用いた結果です．
全データ点を正しく分類できていることがわかります．
その他
SVMのオンライン版というとPassive Aggressiveがありますが，この論文の第1章でPAとの違いが述べられています．
あと，Pegasosって名前めっちゃかっこいい．
(2015/7/8) 予測の式が微妙に間違っていたので修正．
(2015/7/10) 日本語を修正．
(2016/1/15) 引用論文のリンクを修正．