覚え書きブログ

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

Windowspythonで、サポートベクトルマシン(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を実装している。

制約ありのマージン最大化問題を、ラグランジュ双対問題に書き換える。
1)ラグランジュ関数[\tex L(\textbf{w},b,\textbf{a})]を導入する。ここで、[\textbf \textbf{a}]は、[\tex N]次元のラグランジュ乗数である。
2)[\tex L(\textbf{w},b,\textbf{a})]を、[\tex \textbf{w}]と[\tex b]について偏微分してゼロとおき、式を2つ作る。
3)2)の式を、[\tex L(\textbf{w},b,\textbf{a})]に代入し[\tex \textbf{w}]と[\tex b]を消し、停留点(目的関数の負の勾配と制約関数の勾配が釣り合う点)を求める。
4)ラグランジュ関数を最小化するラグランジュ乗数aを求めることにより、重要な制約式(サポートベクトル)を選択する。
f:id:hirotaka_hachiya:20170619194026p:plain
f:id:hirotaka_hachiya:20170619194232p:plain
f:id:hirotaka_hachiya:20170619194319p:plain

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()
    #------------------------------------  

実行結果は、以下のような感じになる。
f:id:hirotaka_hachiya:20170619194434p:plain

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

パワーポイントのアニメーション分割

講義用に資料を配布するときに困るのが、アニメーションで図やテキストが重なって状態で印刷されて見にくくなることだ。
そんなときは、以下のサイトで紹介されているanime2Flipというパワーポイントのアドインが便利。
簡単に指定したスライドのアニメーションを分割してくれる。
http://product-blog.logosware.com/archives/14700

Faster R-CNNの拡張

続々と出てくるFaster R-CNNの改良の論文のなかから、面白いものを挙げていく。

  • Mask R-CNN, Kaiming He et., al arxiv2017

https://arxiv.org/pdf/1703.06870.pdf

  • Subcategory-aware Convolutional Neural Networks for Object Proposals and Detection, Yu Xiang et al., arxiv2017

https://arxiv.org/pdf/1604.04693.pdf