$$ カーネル関数とは、x_nとx_iの近さを表す関数 $$
$$k(\boldsymbol{x},\boldsymbol{x'}) = φ(\boldsymbol{x})^{\mathrm{T}}φ(\boldsymbol{x'}) $$
$$ = φ_1(\boldsymbol{x'})φ_1(\boldsymbol{x'})+φ_2(\boldsymbol{x'})φ_2(\boldsymbol{x'})+ ... +φ_n(\boldsymbol{x'})φ_n(\boldsymbol{x'}) で表される $$
$$ 特徴空間とは、入力ベクトル\boldsymbol{x}の写像 $$
１、特徴空間への写像を考え、その内積をカーネル関数とする
$$k(\boldsymbol{x},\boldsymbol{x'}) = φ(\boldsymbol{x})^{\mathrm{T}}φ(\boldsymbol{x'}) $$
２、カーネル関数を直接定義する
$$k(\boldsymbol{x},\boldsymbol{z}) = (\boldsymbol{x}^{\mathrm{T}}\boldsymbol{z})^2 $$
$$二次元入力(x_1,x_2)を考え、展開すると$$
\begin{align} k(\boldsymbol{x},\boldsymbol{x'}) &= (x_1z_1+ x_2z_2)^2 \\ &= x_1^2z_1^2 + 2x_1z_1x_2z_2 + x_2^2z_2^2\\ &= (x_1^2,\sqrt{2}x_1x_2,x_2^2)(z_1^2,\sqrt{2}z_1z_2,z_2^2)\\ &= φ(\boldsymbol{x})^{\mathrm{T}}φ(\boldsymbol{z}) \end{align}
これは特徴ベクトルの積で表されたので、有効なカーネル関数 $$ちなみに、有効になるための必要十分条件は任意の{x_n}に対して要素がk(x_n,x_n)で与えられるグラム行列kが半正定値であること$$
\begin{align} e^{-γ|x_i-x_j|^2}\\ &=e^{-γ(x_i-x_j)^2}\\ &=e^{-γx_i^2+2γx_ix_j-γx_j^2}\\ &=e^{-γx_i^2-γx_j^2}e^{2γx_ix_j}\\ &=e^{-γx_i^2-γx_j^2}(1+\frac{2γx_ix_j}{1!}+\frac{(2γx_ix_j)^2}{2!}+...)\\ &=e^{-γx_i^2-γx_j^2}(1+\sqrt{\frac{2γ}{1!}}x_i\sqrt{\frac{2γ}{1!}x_j}+...)\\ &= φ(\boldsymbol{x_i})^{\mathrm{T}}φ(\boldsymbol{x_j}) \end{align}
なぜカーネル関数を使うとうまくいくか φという特徴空間への写像を自由に選べるから
例えばgaussian kernel だと、特徴空間の次元が無限になり、超平面（線形関数）で分割できるようになる。 また、カーネル関数を用いると、φがわからずとも内積を計算できる。
この前は、「線形分離がおおよそできる」という条件を満たす任意のφについて成り立つ式変形をしていた。
マージンの外側の領域を0、正しく分類されているがマージン内のものを0〜1、誤分類を1とする、スラック変数を定義する。 $$C\sum_{k=1}^{n} スラック変数 + \frac{1}{2}|\boldsymbol{w}|^2 \\ $$ Cの値を無限大にすれば、ハードマージンとなる
# 非線形SVM（ソフトマージン）
# cvxoptのQuadratic Programmingを解く関数を使用
import numpy as np
from scipy.linalg import norm
import cvxopt
import cvxopt.solvers
from pylab import *
N = 200 # データ数
C = 0.5 # スラック変数を用いて表されるペナルティとマージンのトレードオフパラメータ
SIGMA = 0.3 # ガウスカーネルのパラメータ
# ガウスカーネル
def gaussian_kernel(x, y):
return np.exp(-norm(x-y)**2 / (2 * (SIGMA ** 2)))
# どちらのカーネル関数を使うかここで指定
kernel = gaussian_kernel
def f(x, a, t, X, b):
sum = 0.0
for n in range(N):
sum += a[n] * t[n] * kernel(x, X[n])
return sum + b
if __name__ == "__main__":
# 訓練データをロード
data = np.genfromtxt("classification.txt")
X = data[:,0:2]
t = data[:,2] * 2 - 1.0 # 教師信号を-1 or 1に変換
# ラグランジュ乗数を二次計画法（Quadratic Programming）で求める
K = np.zeros((N, N))
for i in range(N):
for j in range(N):
K[i, j] = t[i] * t[j] * kernel(X[i], X[j])
Q = cvxopt.matrix(K)
p = cvxopt.matrix(-np.ones(N))
temp1 = np.diag([-1.0]*N)
temp2 = np.identity(N)
G = cvxopt.matrix(np.vstack((temp1, temp2)))
temp1 = np.zeros(N)
temp2 = np.ones(N) * C
h = cvxopt.matrix(np.hstack((temp1, temp2)))
A = cvxopt.matrix(t, (1,N))
b = cvxopt.matrix(0.0)
sol = cvxopt.solvers.qp(Q, p, G, h, A, b)
a = array(sol['x']).reshape(N)
print (a)
# サポートベクトルのインデックスを抽出
S = []
M = []
for n in range(len(a)):
if 0 < a[n]:
S.append(n)
if 0 < a[n] < C:
M.append(n)
# bを計算
sum = 0
for n in M:
temp = 0
for m in S:
temp += a[m] * t[m] * kernel(X[n], X[m])
sum += (t[n] - temp)
b = sum / len(M)
print (S, b)
# 訓練データを描画
for n in range(N):
if t[n] > 0:
scatter(X[n,0], X[n,1], c='b', marker='o')
else:
scatter(X[n,0], X[n,1], c='r', marker='o')
# サポートベクトルを描画
# for n in S:
# scatter(X[n,0], X[n,1], s=80, c='c', marker='o')
# 識別境界を描画
X1, X2 = meshgrid(linspace(-2,2,50), linspace(-2,2,50))
w, h = X1.shape
X1.resize(X1.size)
X2.resize(X2.size)
Z = array([f(array([x1, x2]), a, t, X, b) for (x1, x2) in zip(X1, X2)])
X1.resize((w, h))
X2.resize((w, h))
Z.resize((w, h))
CS = contour(X1, X2, Z, [0.0], colors='k', linewidths=1, origin='lower')
for n in S:
print (f(X[n], a, t, X, b))
xlim(-2, 2)
ylim(-2, 2)
show()
# 台湾大学資料 https://www.csie.ntu.edu.tw/~cjlin/talks/kuleuven_svm.pdf
http://aidiary.hatenablog.com/entry/20100503/1272889097