WIndowsへのcvxoptのインストールとSVM

Windowsのpythonで、サポートベクトルマシン(SVM)の二次計画問題を解くのにcvxoptを使ってみた。

以下から、64bitのPython 2.7 versionのAnacondaをダウンロードして、インストール。
https://www.continuum.io/downloads
※ちなみに、python 3.6 versionでは、cvxoptをインストールできなかった。

次に、以下のコマンドでcvxoptをインストール。

> conda install -c omnia cvxopt
> conda install -c omnia scs
> pip install cvxpy

cvxopt + python2.7の準備ができたので、次は、SVMを実装する。
今回は導入編ということで、線形ハードマージンSVMを実装している。

制約ありのマージン最大化問題を、ラグランジュ双対問題に書き換える。
１）ラグランジュ関数を導入する。ここで、は、次元のラグランジュ乗数である。
２）を、とについて偏微分してゼロとおき、式を２つ作る。
３）２）の式を、に代入しとを消し、停留点（目的関数の負の勾配と制約関数の勾配が釣り合う点）を求める。
４）ラグランジュ関数を最小化するラグランジュ乗数aを求めることにより、重要な制約式（サポートベクトル）を選択する。

svmの実装は、下記のブログを参考にした。ただし、今回は線形ハードマージンSVMである。
aidiary.hatenablog.com

#coding:utf-8

# cvxoptを用いたハードマージン線形SVM

import numpy as np
from scipy.linalg import norm
import cvxopt
import cvxopt.solvers
from pylab import *

N = 120  # データ数

# 直線関数
def f(x1, w, b):
    return - (w[0] / w[1]) * x1 - (b / w[1])
    
# 線形カーネル
def kernel_linear(x, y):
    return np.dot(x, y) 
    
kernel = kernel_linear        

if __name__ == "__main__":

    #------------------------------------
    # 学習データを作成
    
    # 各カテゴリのガウス分布を作成し、データを生成
    cls1 = []
    cls2 = []
	
    mean1 = [1, 2]
    mean2 = [-2, -1]
    mean3 = [5, 6]
    cov = [[1,-0.8], [-0.8, 1]]
    
    #cls1.extend(np.random.multivariate_normal(mean1, cov, N*2/6))
    #cls1.extend(np.random.multivariate_normal(mean3, cov, N*1/6))
    cls1.extend(np.random.multivariate_normal(mean1, cov, N/2))
    cls2.extend(np.random.multivariate_normal(mean2, cov, N/2))
        
    # データ行列Xを作成
    X = vstack((cls1, cls2))
    
    # ラベルベクトルyを作成
    y = []
    for i in range(N/2):
        y.append(1.0)   # クラス1
    for i in range(N/2):
        y.append(-1.0)  # クラス2
    y = array(y)
    #------------------------------------
	
    #------------------------------------    
    # ハードマージンSVMの双対問題を凸二次計画法で解き、ラグランジュ未定乗数を求める
    #
    # 【ハードマージンSVMの双対問題】
    #   min 1/2 a^\top (y^\top x^\top x y) a - 1\top a
    #   s.t. -a \ge 0
    #        a^\top y = 0
    #  を下記の凸二次計画法の各変数（P, q, G, h, A, b）に当てはめる。
    #
    # 【凸二次計画法】
    #   min 1/2 X^\top P X + q^\top X
    #   s.t. Gx \ge h
    #        Ax = b
    #
    #  cvxopt.solvers.qp(P,q,G,h,A,b)

    # Pの作成
    Ptmp = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            Ptmp[i, j] = y[i] * y[j] * np.dot(X[i], X[j])           
    P = cvxopt.matrix(Ptmp)
    
    # q,G,h,A,bを作成
    q = cvxopt.matrix(-np.ones(N))             # 全ての要素が-1の列ベクトル（Nx1）
    G = cvxopt.matrix(np.diag([-1.0]*N))       # 対角成分が-1行列（NxN）
    h = cvxopt.matrix(np.zeros(N))             # 全ての要素が0の列ベクトル（Nx1）
    A = cvxopt.matrix(y, (1,N))                # 各要素がカテゴリラベルの行ベクトル（1xN）
    b = cvxopt.matrix(0.0)                     # スカラー0.0
    
    # 凸二次計画法を解いてラグランジュ未定乗数aを求める
    sol = cvxopt.solvers.qp(P, q, G, h, A, b)
    a = array(sol['x']).reshape(N)             # 'x'がaに対応する
    #------------------------------------

    print("a:",a)

    #------------------------------------    
    # 分類境界w, bをサポートベクトルから求める

    # サポートベクトルのインデックスを抽出
    S = []
    for i in range(len(a)):
        if a[i] < 0.00001: continue
        S.append(i)
    
    # wを計算
    w = np.zeros(2)
    for n in S:
        w += a[n] * y[n] * X[n]
    
    # bを計算
    sum = 0
    for n in S:
        temp = 0
        for m in S:
            temp += a[m] * y[m] * kernel(X[m], X[n])
        sum += (y[n] - temp)
    b = sum / len(S)
    
    print("S:",S)
    print("b:",b)
    #------------------------------------    
    
    #------------------------------------    
    # 学習データを描画

    # カテゴリ1
    x1, x2 = np.array(cls1).transpose()
    plot(x1, x2, 'o', color="#FFA500", markeredgecolor='k', markersize=14)

    # カテゴリ2    
    x1, x2 = np.array(cls2).transpose()
    plot(x1, x2, 's', color="#FFFF00", markeredgecolor='k', markersize=14)
    
    # サポートベクトルを描画
    for n in S:
        plot(X[n,0], X[n,1], 'o', color='none', markeredgecolor='r', markersize=18, markeredgewidth=3)
    
    # 識別境界を描画
    x1 = np.linspace(-6, 6, 1000)
    x2 = [f(x, w, b) for x in x1]
    plot(x1, x2, 'g-')
    
    xlabel("x1", fontsize=14)
    ylabel("x2", fontsize=14)
    tick_params(labelsize=14)

    xlim(-6, 6)
    ylim(-6, 6)
    show()
    #------------------------------------  

実行結果は、以下のような感じになる。

次回は、カーネルモデルに拡張する。