lazy diary

統計とその周辺

ハードマージンSVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回はハードマージンSVMというものを実装します。
ソフトでもない上に非線形でもないという最も基本的なSVMです。ハードマージンSVMさえ理解してしまえば他のSVMはその拡張みたいなものです。

必要なものをimportしておきます。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

このcvxoptというのは二次計画問題のソルバーです。詳しくはこの記事の最後にあるブログを参考にしてください。google colabにもデフォルトではインストールされていないようなので、

pip install cvxopt

で導入しておきましょう。

データは相変わらずいい加減に作ったものです。とは言いつつも今回はハードマージンSVMですので、線形分離できないデータは困ります。かならず超平面で分離できるデータをチョイスしてください。2次元なら直線で完全に分離できるものです。

#create data
x_1 = 3
y_1 = 7
x_2 = 7
y_2 = 2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

ここからが本編です。流れとしては
1. データから二次計画問題を解く
2. パラメータwを求める
3. サポートベクトル集合を求める
4. サポートベクトル集合から切片bを求める
となっています。

class HardMarginSVM():
  
  def __init__(self, X, y):
    self.X = X
    self.y = y
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    G = matrix((-1) * np.identity(num_of_data))
    h = matrix(np.zeros((num_of_data, 1)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon:
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

これを上のデータに当てはめてみましょう。

svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022105408p:plain
黄色になっている部分はサポートベクトルです。SVMの特徴として実質的にサポートベクトルにしか依存しないので、次のような外れ値が入っても影響がありません。

data.loc[100] = [0, 20, 1]
svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022110158p:plain

参考

今回、二次計画問題を解くにあたって参考にしました。SMOアルゴリズムについてはまだ理解できていないため、今回はcvxoptというソルバーに頼りました。
二次計画法(Quadratic Programming)の概要とPythonソルバーcvxoptの使い方 - MyEnigma