14. AI技術:Support Vector Machine 入門

異なる色の領域を隔てる直線的な構図の写真
Photo by Craig Whitehead on Unsplash

Support Vector Machine(SVM)は,データを分ける境界を「最も余裕のある位置」に置く分類手法である.この回では,境界・マージン・support vectorの直観から始め,Juliaで分類器を作り,最後にカーネルを用いた非線形分類まで見通す.

この回の到達目標

  1. 境界を理解する超平面で二種類を分ける考え方を説明できる
  2. マージンを理解するsupport vectorが境界を決める理由をつかむ
  3. 分類器を作るOptimとLIBSVMを使った例を実行できる
  4. 非線形へ拡張するsoft marginとkernelの役割を説明できる

1. Support Vector Machine(SVM)とは

Support Vector Machine(以降 SVM)とは,基本的には,適切に学習することでデータ点を「二種類1に分類できるようになる」ことが目的の分類器(classifier)である AI の一種で,以下のような仕組みになっていると考えれば良い.

SVMのおおまかなアイディア

1二種類の既知データ点を空間に配置した図既知データを配置特徴量を座標として,二種類の点を空間に置く
2二種類のデータ点の間に分類境界を置いた図境界を学習両側の余裕が大きくなる超平面を求める
3未知のデータ点を境界の側によって分類する図未知データを分類境界のどちら側にあるかで種類を判定する
図:SVMの基本は「配置する → 境界を学ぶ → 未知データを判定する」の3段階である.
  1. 空間上に既知のデータ点を配置する. このとき,データ点がおおまかに2つに別れて配置されている,と想定する.
  2. 一つの超平面を用意し「その両側で空間を二分して」,そのどちらにデータ点があるかでデータを 2つの種類にうまく分類できるようにする.言い換えると,そのような適切な超平面を既知のデータ点から計算する.これが SVM の学習過程である.このとき,この超平面に一番近いデータ点を Support Vector 2と呼ぶ.
  3. 未知のデータ点に対し,上で求めた超平面のどちらにあるかで,分類を判断する.

先に押さえる二つの疑問

さて,アイディアはシンプルで分かりやすい. ただ,アイディアを理解した時点で根本的な疑念がいくつか生じるだろう. そのあたりを先に記述しておこう.

  1. 疑念.     空間にデータ点を配置しても,超平面でうまく二分できるようにデータ点が配置されているとは限らないのでは?

    回答.     そのとおり.そのため,通常はデータ点を

       ・非線形変換する
       ・(上の変換後のデータと組み合わせて)データ点の配置空間の次元を上げる

    ことでうまく二分できる「可能性を上げる」. これについては, SVM with polynomial kernel visualization (HD) (Youtube 動画, by udiprod)が大変わかりやすいイメージ動画になっているので一度見ておこう.
    さらに,

       ・うまく分離できないデータ点の存在を許すソフトマージンを採用する3

    ことで全体の柔軟性を上げるのだ. より詳しいことは後述する.

  2. 疑念.     「適切な」超平面を計算する手順は? とくに上記の非線形変換などを行うと,うまくできる気がしないが…

    回答.     「なるべく良い」超平面を探す最適化問題に帰着させることとちょっとした工夫で,既存の知識で可能になる. 重要なポイントを書いておこう.

       ・分類失敗するデータ点も許す「ソフトマージン」で超平面を計算する!

       ・非線形変換を行う場合,最適化問題の主問題ではなく双対問題を解く

       ・この双対問題はデータ点の二次形式4になる

       ・次元やデータ点数が多くなると二次形式の計算量増大がつらい

       ・二次形式相当の結果をいきなり計算できるカーネルを採用することで,計算量増大をゆるやかにできる.これはカーネルトリックと呼ばれる重要な工夫だ.

詳しくは後述する.

margin境界からの余裕大きいほど未知データに強い境界を期待できる
soft margin少数の分類ミスを許す現実の重なりや外れ値に対応する
kernel非線形な境界へ高次元での内積を直接計算する
図:SVMを実用的な分類器にする三つの考え方.

2. 線形な境界でほぼ分離できる場合

データ点を素直に空間に配置しただけで二種類のグループに超平面で分離が可能な場合から考えよう. こんな都合の良い場合は少ないだろうが,アイディアの根本なのでここから始めよう. ちなみにこれを「線形分離可能(linearly separable)な場合」と呼ぶ.

ハードマージン:全サンプルをきれいに分離する

この場合,以下のような図で全体を考えることができる.

分類境界,二本のmargin境界,support vector,法線ベクトルwの関係を示す図
図:境界に最も近い点がsupport vectorであり,その点までの距離がmarginを決める.

ただし,$\boldsymbol{w}, b$ はそれぞれこの超平面の法線ベクトルと切片で,超平面自身(上の図中の水色の線)が

\begin{equation} \boldsymbol{w}\cdot\boldsymbol{x} - b = 0 \end{equation}

という方程式で書かれるとしている.

そして上の図で左上側にあるグループ「青」を $X_{+}$, 右下側グループ「赤」を $X_{-}$ と書くことにして, それぞれ相手グループ側に一番近い点 $\boldsymbol{x}_{\pm \mbox{c}}$ として

\begin{equation} \left\{ \begin{array}{rcl} \boldsymbol{w}\cdot\boldsymbol{x_{+\mbox{c}}} - b & = & 1, \cr \boldsymbol{w}\cdot\boldsymbol{x_{-\mbox{c}}} - b & = & -1 \end{array} \right. \label{eq:on-margin-line} \end{equation}

を満たすとしよう.

  超平面の情報,すなわち $\boldsymbol{w}$, $b$ を決めて計算すると,ぎりぎりに乗っている $\boldsymbol{x_{\pm\mbox{c}}}$ に対して $\boldsymbol{w}\cdot\boldsymbol{x_{\pm\mbox{c}}} - b$ を計算しても 1 や -1 ちょうどにならないのでは? と思うだろう. 実はそのとおりだ.
実は $\boldsymbol{w}$, $b$ には定数倍の自由度があり,それが原因だ. 具体的には,非ゼロ定数 $C$ を用いて $C\boldsymbol{w}$, $Cb$ という法線ベクトルと切片を使った超平面: $(C\boldsymbol{w})\cdot\boldsymbol{x} - (C b) = 0$ は,元の超平面: $\boldsymbol{w}\cdot\boldsymbol{x} - b = 0$ と全く 同じ超平面を指すのだ. だから,上の話は「この定数 $C$ を調整して $\boldsymbol{w}\cdot\boldsymbol{x_{\pm\mbox{c}}} - b$ がちょうど $\pm 1$ になるように $\boldsymbol{w}$, $b$ を定数倍する操作が隠されている」 という意味だ,と思って良い.
ちなみに望ましい $C$ の計算方法だが,これは簡単で,今の $\boldsymbol{w}$, $b$ を用いて $1/C := \boldsymbol{w}\cdot\boldsymbol{x_{\pm\mbox{c}}} - b$ と求めれば良い.そして超平面のパラメータの正しい値を $C\boldsymbol{w}$, $Cb$ にすれば OK だ.

さて,話を戻そう. グループの他の点も考えると,上の条件は

\begin{equation} \left\{ \begin{array}{rcll} \boldsymbol{w}\cdot\boldsymbol{x} - b & \geq & 1 & \mbox{ for all } \boldsymbol{x} \in X_{+}, \cr \boldsymbol{w}\cdot\boldsymbol{x} - b & \leq & -1 & \mbox{ for all } \boldsymbol{x} \in X_{-} \end{array} \right. \label{eq:in-group} \end{equation}

と書き直せる.

そして,このグループ間のマージン(図を見よ)は

\begin{equation} \frac{2}{\| \boldsymbol{w} \|} \end{equation}

となるので,これをハードマージンと呼び,分類がなるべくしやすいように,この値をなるべく大きくする超平面を求めることにする,というのがストーリーだ.

つまり, 式(\ref{eq:on-margin-line})と式(\ref{eq:in-group})を満たす範囲でなるべく $\| \boldsymbol{w} \|$ を小さくすればよいという,最小化問題なのだ.

さて,この問題をもう少し整理しよう. まず,各データ点 $\{\boldsymbol{x}_i\}_i^N$ の分類を値 $\{y_i\}_i^N$で表すことにし,例えば左上のグループ「青」はすべて $y_i = +1$, 右下「赤」は $y_i = -1$ であるとしよう.

すると式(\ref{eq:on-margin-line})と式($\ref{eq:in-group}$) はあわせて \begin{equation} y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 1 \,\,\,\mbox{ for all } (\boldsymbol{x}_i, y_i), \label{eq:cond-hardmargin} \end{equation}      ただし,上の式の等号を満たす $y_i = 1, y_j = -1$ が存在する.

という条件に書き直せるので,ハードマージン下での問題全体は

  制約条件付き最小化問題:
      $\| \boldsymbol{w} \|$ を最小化せよ.
      ただし,全てのデータ $(\boldsymbol{x}_i, y_i)$ に対して $y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 1$ を満たし, かつ,この式の等号を満たす $y_i = 1, y_j = -1$ が存在するようにせよ.

を解けば良い,ということになる. これは典型的な最小化問題なので,あとは最適化問題の技術や知見,ライブラリに頼れば良いだけだ.

なお,上の最小化問題を解くと超平面の法線ベクトル $\boldsymbol{w}$ が求まる. ならば超平面の切片 $b$ は? と思うかもしれないが,それは $\boldsymbol{w}$ と, $\boldsymbol{x}_{\pm\mbox{c}}$ のどれかが特定できていれば 式($\ref{eq:on-margin-line}$)から計算できるので問題ない.

学習後の分類

上の最小化問題を解いて $\boldsymbol{w}, b$ が求まっていれば,未知のデータ $\boldsymbol{x}$ の分類は簡単だ.

$y = \boldsymbol{w}\cdot\boldsymbol{x} - b$ を計算して, $y$ がプラスならばグループ「青」$X_{+}$ に分類し, $y$ がマイナスならばグループ「赤」$X_{-}$ に分類すれば OK! だ.



ソフトマージン:多少のミスを許す

そのまま配置で「ほぼ」分離できるケースに問題を緩和して考えよう.

つまり,全てのデータ点に対して $y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 1$ を要求していたハードマージンに対し, 多少は条件を破っても良いとするのだ.

これがソフトマージンである.

全体には,

  • 条件を満たしていないデータ点について「ペナルティ量」を用意し,
  • 本来小さくしたい量 $\| \boldsymbol{w} \|$ にペナルティを足して小さくすべき全体量を定義する

という仕組みを考えるのだ.

marginを広くするλ ||w||2境界に余裕を持たせる
違反を少なくするhinge loss誤分類やmargin内の点にペナルティ
両者の折り合いsoft margin現実のデータに合う境界を探す
図:soft marginは,marginの広さと分類違反の少なさを一つの目的関数で調整する.

より詳しく解説しよう.

まず,本来守るべき制約条件 $y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 1$ $\Leftrightarrow$ $0 \geq 1 - y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b)$ に対し, この条件を守れない場合のペナルティを,

\begin{equation} 1 - y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \end{equation}

が 0以上になった場合その数値にする,というアイディアである(条件を守れていない場合はこの値がプラスになる).

そして,条件を守れている場合はもちろんペナルティはゼロで良いので,両方の場合を合わせて,ペナルティを

\begin{equation} \mbox{ペナルティ項 } := \max( 0, 1 - y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) ) \end{equation}

とすれば扱いやすい. というわけで,ペナルティをこう設定しよう.

あとは,本来小さくしたい量である $\| \boldsymbol{w} \|$ と合わせて考えれば良いので,重み $\lambda$ をつけることにして,たとえば次のような量 $L$

\begin{equation} L := \lambda \| \boldsymbol{w} \|^2 + \frac{1}{N} \sum_{i=1}^N \left\{ \max( 0, 1 - y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) ) \right\} \end{equation}

を考えて, 式(\ref{eq:on-margin-line})を満たす点 $\boldsymbol{x_{+\mbox{c}}}$, $\boldsymbol{x_{-\mbox{c}}}$ が存在するように $\boldsymbol{w}$ と $b$ を調整したうえで, これを最小化すればいい,ということになる. つまり,

  最小化問題:
      与えられた正係数 $\lambda$, データ $(\boldsymbol{x}_i, y_i)_i^N$ に対し, 式($\ref{eq:on-margin-line}$) を満たす点 $\boldsymbol{x_{+\mbox{c}}}$, $\boldsymbol{x_{-\mbox{c}}}$ が存在するように $\boldsymbol{w}$ と $b$ を調整したうえで, 量 $L$ を最小化する $\boldsymbol{w}, b$ を求めよ.

を解いて $\boldsymbol{w}, b$ を求めれば良いのだ.

ちなみに,重み係数 $\lambda$ を大きくすると本来の小さくしたい目的量 $\| \boldsymbol{w} \|^2$ を重視し,ペナルティが軽視される. 逆に係数 $\lambda$ を小さくするとペナルティが重視され,マージンをハードマージンにとった状況に近くなる.

3. 実習1:Optimで線形SVMを作る

さて,このあたりで簡単な例で実際にやってみよう. 分離が可能な簡単な例で,ソフトマージンで超平面を計算してみるのだ(ハードマージンでも良いぞ).

なお,最適化のために Optim package を用いよう.

Optimパッケージのインストール(すでにインストールされている場合は不要)

おそらく多くの人は未インストールだろうから、Optim package を下記のようにしてインストールしておこう.

1using Pkg
2Pkg.add("Optim")

データを用意する

さて,話を戻して,まずは準備だ.

1using Optim
2using LinearAlgebra
3using Plots

次に,対象データを作成しよう. まあ何でも良いのだが,今回は下の図を考える.

1a, b = 3.0, 3.5
2fu(x) = (x-a)^2 + b
3fd(x) = -(x-b)^2 + a
4
5X = 0:0.01:5
6plot(X, fu)
7plot!(X, fd)
青い上側曲線と赤い下側曲線でデータ生成領域を示したグラフ
図:青い曲線より上を正例,赤い曲線より下を負例の生成領域とする.

そして,この青い線より上に「青」のデータを,赤い線より下に「赤」のデータを用意する.ちょっと長いが次のようにすればよいだろう. まず,線の上下にランダムにデータを取る関数を次のように作成する.

 1L = 20.0
 2function upper()
 3    x =  L * rand() / 4
 4    y = 0.0
 5    while y < fu(x)
 6        y = L * rand()
 7    end
 8    return [x, y]
 9end
10
11function downer()
12    x = L * rand() / 4
13    y = 10.0
14    while y > fd(x)
15        y =  L * ( rand() - 0.5 )
16    end
17    return [x, y]
18end

そして,2000点ほどデータ点を作ろう. $X_p$ が青のデータ点集合(正値集合), $X_n$ が赤のデータ点集合(負値集合)だ.

1N = 2000
2
3Xp = upper()
4Xn = downer()
5
6for i in 1:div(N,2)-1
7    Xp = hcat( Xp, upper() )
8    Xn = hcat( Xn, downer() )
9end

どんな点配置か,プロットして確認しよう.

1scatter(Xp[1,:], Xp[2,:])
2scatter!(Xn[1,:], Xn[2,:])
青と赤の二群に分かれた2000個のサンプル点
図:生成した正例と負例.この例は直線でほぼ分離できる.

次に,なにも調整していない $\boldsymbol{w}, b$ を正しく調整する函数を用意しよう. なにをどうするかのイメージは下図のとおりだ.

法線ベクトルと切片を定数倍してmargin境界を調整する考え方の図
図:同じ超平面を表す定数倍の自由度を使い,最も近い点で値が1になるよう調整する.

そのために,まずは \begin{equation*} y_i (\boldsymbol{w}\cdot\boldsymbol{x_i} - b) \end{equation*} の値を全て調べて,その中で「正な値のなかで最小」な値(上の図の c )を求める函数を下記のように用意する.

 1function r_plusmin(w, b)
 2  min_r = 100000.0
 3  len = size(Xp)[2]
 4
 5  for i in 1:len
 6    r =  w  Xp[:,i] - b
 7    if (r > 0.0) && (min_r > r)
 8      min_r = r
 9    end
10    r = - ( w  Xn[:,i] - b )
11    if (r > 0.0) && (min_r > r)
12      min_r = r
13    end
14  end
15
16  return min_r
17end

  上で r = の式の右辺で wXp などの間にある記号 ⋅(見えにくいが)は,ベクトルの内積を意味する記号で,\cdot & tabキー で入力できる. dot 関数の syntax sugar だ.

この結果を用いて,まずはギリギリの点で上の式が 1 になるように $\boldsymbol{w}, b$ を調整するのだ. 次のような函数を用意すれば良い.

1function tune(w, b)
2  c = r_plusmin(w, b)
3  return w ./ c,  b/c
4end

この出力を $\tilde{\boldsymbol{w}}, \tilde{b}$とすると,これはギリギリの点での条件を満たす,いわば調整済みのパラメータになる. この値を使ってペナルティ等を計算すれば良い.

そして,あと必要なのはペナルティを計算して,$\boldsymbol{w}$ の大きさと足し合わせて量 $L$ を作ることだ. 次のように定義すればよいだろう.

 1function func_L(v)
 2  λ = 0.01
 3  w = v[1:2]
 4  b = v[3]
 5  
 6  # まず最初にきちんと調整する
 7  w_true, b_true = tune(w,b) 
 8
 9  # penalty の計算
10  len = size(Xp)[2]
11  penalty = 0.0
12
13  for i in 1:len
14    r =  w_true  Xp[:,i] - b_true
15    if r < 0.0
16      penalty += 1 - r
17    end
18    r = - ( w_true  Xn[:,i] - b_true )
19    if r < 0.0
20      penalty += 1 - r
21    end
22  end
23  
24  return λ * (norm(w_true)^2) + penalty / (2*len)
25end

あとは $\boldsymbol{w}, b$ の値を変えて $L$ が一番小さくなるところを探すのだ. この「最小点を探す作業」は Optim package におまかせでいい.次のような感じだ. ちなみに探索初期値は,超平面が右上がりになっている,というぐらいの設定だ.

1res = optimize(func_L, [-1.0, 1.0, 1.0 ])
2# 小さくしたい関数と,探索初期値を渡す.
 * Status: success

 * Candidate solution
    Final objective value:     3.590264e-02

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    163
    f(x) calls:    322

どうやら成功したようだ.結果は以下のようにして取り出せる.

1v = res.minimizer
2w, b = v[1:2], v[3]
([-0.15465492690086702, 0.8442553526118208], 2.205029125079519)

ただしこの $\boldsymbol{w}, b$ は調整してない値だ.調整しておこう. 皆の結果も似たような数字になるはずだ(データを乱数で作っているので,やるたびに少し違う)

1w_true, b_true = tune(w,b)
([-0.34141763088297605, 1.8637858367988631], 4.867842460623156)

分離具合をプロットしてみてみよう.

1scatter(Xp[1,:], Xp[2,:])
2scatter!(Xn[1,:], Xn[2,:])
3
4plot!( x -> ( -w_true[1]*x + b_true )/w_true[2], label = "hyperplane" )
青と赤のデータ点を学習済みの直線で分離した結果
図:Optimで求めた超平面が二群を分離している.

ふむ,たいへんうまく分離できているようだね. これで OK だ.

ちなみに,作った超平面に基づく分類器のマネごとも作ってみよう. 次のような感じだね.

1classifier(x) = sign( w_true⋅x - b_true)

試しに,「青」になりそうなデータ値を判定させてみる.

1x_un = [0.5, 5.0]
2classifier(x_un)
1.0

正の値を返したので,たしかに「青」と分類しているね.

同様に「赤」っぽいデータ値も判定させてみよう.

1x_un = [0.5, -5.0]
2classifier(x_un)
-1.0

これも負値を返したので正しく「赤」と分類している.

というわけで,SVM による分類器が正しく作れた,ということになるね. 簡単な問題の場合はこんな感じだ.

4. 線形な境界では分離が難しい場合

さて,これから示す問題がほんとうにやりたい問題だ. データ点を配置しても素直に超平面で分離なんかできないぞ,という,よくありそうなケースを考える.

なお,以降の話を含めて2分で解説している 動画@youtube があるので見るのも良いだろう.

さてこの場合,おおよそ主流の組み合わせが決まっていて,

      非線形変換 & 双対問題を考える & カーネルトリックの採用

という形になっている5. ここでもこの流れを紹介しよう.

元の空間x直線では分けにくい
非線形変換φ(x)分けやすい高次元空間を考える
内積だけ必要φ(xi)・φ(xj)変換後の座標そのものは不要
kernelで近道k(xi, xj)内積の結果を直接計算する
図:カーネルトリックでは,高次元の座標を明示的に作らず,必要な内積だけをkernelで求める.

非線形変換

まず非線形変換であるが,これは形式的には \begin{equation} \boldsymbol{x} \mapsto \boldsymbol{\varphi}(\boldsymbol{x}) \end{equation} とするだけだ. ただし,このあとカーネルトリックを採用するため,
この$\boldsymbol{\varphi}$ の形は陽に与えないでよい6

カーネルトリック

非線形変換後の双対問題では,$\boldsymbol{\varphi}(\boldsymbol{x})$ が単独ではなく,二つの変換済みベクトルの内積として現れる.そこで,この内積を二変数関数

\begin{equation} k(\boldsymbol{x}_i,\boldsymbol{x}_j) \end{equation}

で直接計算する.これがカーネルトリックである.高次元の座標を実際に作らなくても非線形な分類境界を利用できる点が重要だ.双対問題の式と分類式の導出は,後半の「発展・補足編」にまとめる.

代表的なカーネル

ちなみにカーネルの種類等であるが,まあ典型的なのは次のような感じだ. よくわかんないときはとりあえず G-RBF を使っておけばいいんじゃないかな.

  • カーネルの例
    名前 関数 $k( \boldsymbol{x}_i, \boldsymbol{x}_j ) = $ 備考
    多項式 $( \boldsymbol{x}_i\cdot\boldsymbol{x}_j + r )^d$ 意外と有効だ
    ガウシアン動径基底関数
    (Gaussian RBF, Gaussian radial basis function)
    $\displaystyle \exp\left( - \frac{ \| \boldsymbol{x}_i - \boldsymbol{x}_j \|^2 }{2\sigma^2} \right)$ 広く使われる
    sigmoid $\tanh(\kappa \boldsymbol{x}_i\cdot\boldsymbol{x}_j + c )$ ただし,$c < 0 < \kappa$



5. 実習2:LIBSVMでアヤメを分類する

ここでは SVM を扱うちょっと有名なライブラリ LIBSVM7 を julia で使えるようにした LIBSVM package を採用して実習してみよう.

LIBSVMパッケージのインストール

LIBSVM パッケージを未インストールの場合は下記のようにしてインストールしておこう(VIORAS環境では既にインストール済みだ).

1using Pkg
2Pkg.add("LIBSVM")

  ちなみに LIBSVM package の使い方だが,using LIBSVM してある状況で ?svmtrain とするとデフォルト値がどうなっているかとか使い方が分かるよ. デフォルトはガウシアン RBF をカーネルにした非線形 SVM だね. あとは LIBSVM の情報を見ると良いね.

では,まずは典型的な例として,前回やったアヤメのデータを分類させてみよう. 下記のように準備 & データのダウンロードをして,

 1using LIBSVM
 2using MLDatasets
 3using Statistics
 4using DataFrames
 5
 6iris = Iris( as_df = false )
 7# 利用データをダウンロード.一回だけダウンロードすれば良い.
 8
 9features_raw = iris.features
10labels_raw   = reshape(iris.targets, :, 1)
11# 花の4つの数字情報と,花の種類を配列に入れる.
12# 扱いやすいよう,種類の配列を転置しておく.

今回はデータの半分を SVM の学習用データに,残りの半分を SVM の性能評価に回そう. 次のようにデータを半分に分けておく.

1# データ点
2X_train = features_raw[:, 1:2:end ] # 学習用
3X_check = features_raw[:, 2:2:end ] # 検証用
4
5# 分類.この場合はアヤメの種類.3種類ある.
6y_train = labels_raw[ 1:2:end ] # 学習用
7y_check = labels_raw[ 2:2:end ] # 検証用
全データ150個4特徴・3品種
奇数番目75個学習用:境界を決める
偶数番目75個検証用:未知データとして測る
図:同じデータで学習と評価をしないよう,アヤメ150個を半分ずつに分ける.

そうそう,アヤメの花は 3種類あったので,上の SVM のストーリーからずれる. どうするの? と思うだろうが,2種類判定を繰り返せばよいのだ. そして,以下のように,LIBSVM package はそのあたりを自動でやってくれる.

  LIBSVM package では,分類が 2種類より大きい $n$ 種類になっても大丈夫.2種類分類器を ${}_n C_2$ 個作ってそれぞれで判定し,多数決で判定するタイプの多種類分類(one-against-one strategy と呼ぶ)を自動でやってくれる.

setosavsversicolor
setosavsvirginica
versicolorvsvirginica
3分類器の多数決→ 最終的な1品種
図:3種類なら ${}_3C_2=3$ 個の二値分類器を作るone-against-one方式になる.

で, LIBSVM ではいきなり以下の 1行で SVM を学習,構築できる.

1model = svmtrain( X_train, y_train )
2# SVM の学習.超平面のパラメータが計算される.
LIBSVM.SVM{InlineStrings.String15, LIBSVM.Kernel.KERNEL}(SVC, 
LIBSVM.Kernel.RadialBasis, nothing, 4, 75, 3, 
InlineStrings.String15["Iris-setosa", "Iris-versicolor", "Iris-virginica"], 
…略…
0.0, [0.8070380364343522 0.8047780005392675; 0.0 0.741676866696239; … ; -0.0 -1.0; -0.0 -1.0], 
Float64[], Float64[], [0.12550506116227908, 0.21415319269230185, -0.07664844428517169], 
3, 0.25, 200.0, 0.001, 1.0, 0.5, 0.1, true, false)

  この例では,前回の Deep Learning に比べて SVM の学習が一瞬で終わることに注目しよう.問題規模やモデルによって事情は異なるので一概には比較できないが,小規模な表形式データではSVMが高速かつ高精度に働くことがある.

さて,出来上がった SVM分類器を用いて,検証用データを分類させてみよう.

1y_predict = svmpredict( model, X_check )[1]
2# 検証用データを分類させてみる.
75-element Vector{String}:
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
 "Iris-setosa"
…略…

そして,分類結果を正解と比較してみよう.たくさんあるので,正解率をいきなり計算する形でいいかな.

1mean( y_predict .== y_check )
2# 検証データの分類結果と,正解を比較して正解平均率を求める.
0.9333333333333333

ほう,いきなり 93% の正解率を出したか. データの正規化等を行っていないし,学習に用いたデータはわずか 75個にすぎないのにだいぶ良い値だ. こうした分類に関してはやはり SVM は良い性能を持つのだな,と実感できる.

検証用データ93.3%70 / 75個に相当
図:学習に使っていない75個に対する正解率.

発展・補足編

双対問題をもう少し詳しく見る

そもそも最適化問題の理論では,もとの問題(主問題 primal problem)に対して,双対問題 dual problem と呼ばれる問題を考える.

双対定理(duality theorem)

適切な凸性や制約条件など,強双対性が成り立つ条件のもとでは,主問題と双対問題の最適値が一致する.そのため,両者の「難易度」が異なる場合は,扱いやすい方を解けばよいという考え方につながる.

詳細1は省くが,上のソフトマージンの最小化問題にも双対問題を考えることができ,しかもこちらの方が非線形分類へ拡張しやすい.以下のようになる.

  上の最小化問題の双対問題:
      与えられた正係数 $\lambda$, データ $(\boldsymbol{x}_i, y_i)_i^N$ に対し,

\begin{equation} f(c_1, c_2, \cdots, c_N) = \sum_{i=1}^N c_i - \frac{1}{2}\sum_{i=1}^N \sum_j^N y_i c_i \boldsymbol{x}_i\cdot\boldsymbol{x}_j y_j c_j \label{eq:normal-dual-problem} \end{equation}

      を最大化する $\{c_i\}_i^N$ を求めよ.
      ただし,$\sum_i^N c_i y_i = 0$ と全ての $c_i$ に対して $0 \leq c_i \leq 1/(2N\lambda)$ を満たせ.

そして,超平面の情報であるが,これを解いて求まる $\{c_i\}_i^N$ によって

\begin{equation} \boldsymbol{w} = \sum_{i=1}^N c_i y_i \boldsymbol{x}_i \end{equation}

と法線ベクトルが定まる.切片 $b$ についてはハードマージンのときと同じだ.

非線形変換後の双対問題と分類式

次に双対問題を考える.これは,ソフトマージンの双対問題に非線形変換を組み込んだ,次の問題になる.

  非線形変換時の双対問題(カーネルトリック採用前):
      与えられた正係数 $\lambda$, データ $(\boldsymbol{x}_i, y_i)_i^N$ に対し,

\begin{equation} f(c_1, c_2, \cdots, c_N) = \sum_{i=1}^N c_i - \frac{1}{2}\sum_{i=1}^N \sum_j^N y_i c_i \boldsymbol{\varphi}(\boldsymbol{x}_i)\cdot\boldsymbol{\varphi}(\boldsymbol{x}_j) y_j c_j \label{eq:nonlinear-dual-problem} \end{equation}

      を最大化する $\{c_i\}_i^N$ を求めよ.
      ただし,$\sum_i^N c_i y_i = 0$ と全ての $c_i$ に対して $0 \leq c_i \leq 1/(2N\lambda)$ を満たせ.

ここで,非線形変換関数 $\boldsymbol{\varphi}$ が最適化問題の中で単独では現れず,必ず二次形式ないしはベクトルの内積

\begin{equation} \boldsymbol{\varphi}(\boldsymbol{x}_i)\cdot\boldsymbol{\varphi}(\boldsymbol{x}_j) \label{eq:ip} \end{equation}

の形で現れることに着目する.この項を,カーネル(kernel)と呼ばれる「計算量が少ない二変数関数」

\begin{equation} k( \boldsymbol{x}_i, \boldsymbol{x}_j ) \label{eq:kernel} \end{equation}

で置き換える.これがカーネルトリックである.**再生核Hilbert空間(reproducing kernel Hilbert space)**の議論によって,一定の前提条件を満たせばこの置き換えが可能だと分かっている.こうすることで,高次元での二次形式や内積を明示的に計算せずに済む.

したがって,最終的には次の問題を解くことになる.

  非線形変換 & カーネルトリック採用時の双対問題:
      与えられた正係数 $\lambda$, データ $(\boldsymbol{x}_i, y_i)_i^N$ に対し,

\begin{equation} f(c_1, c_2, \cdots, c_N) = \sum_{i=1}^N c_i - \frac{1}{2}\sum_{i=1}^N \sum_j^N y_i c_i k( \boldsymbol{x}_i, \boldsymbol{x}_j ) y_j c_j \label{eq:nonlinear-kernel-trick-dual-problem} \end{equation}

      を最大化する $\{c_i\}_i^N$ を求めよ.
      ただし,$\sum_i^N c_i y_i = 0$ と全ての $c_i$ に対して $0 \leq c_i \leq 1/(2N\lambda)$ を満たせ.

そして超平面の決定と分類では,超平面の法線ベクトル $\boldsymbol{w}$ を明示的に計算しなくてよい.次の変形から,カーネルだけあればよいことが分かる.

  • 超平面の切片 $b$ の計算:
    マージンぎりぎりにある $\boldsymbol{x}_{\mbox{c}}$ を一つ見つけておいて,次のように計算する.

\begin{eqnarray} b = \boldsymbol{w}\cdot\boldsymbol{\varphi}(\boldsymbol{x}_{\mbox{c}}) - y_{\mbox{c}} & = & \sum_{i=1}^N c_i y_i \boldsymbol{\varphi}(\boldsymbol{x}_i)\cdot\boldsymbol{\varphi}(\boldsymbol{x}_{\mbox{c}}) - y_{\mbox{c}} \nonumber \cr & = & \sum_{i=1}^N c_i y_i k( \boldsymbol{x}_i, \boldsymbol{x}_{\mbox{c}} ) - y_{\mbox{c}} \end{eqnarray}

  • 未知のデータ $\boldsymbol{x}$ の分類:
    次の量がプラスならば $X_{+}$ の元,マイナスならば $X_{-}$ の元であると判定する.

\begin{equation} \boldsymbol{w}\cdot\boldsymbol{\varphi}(\boldsymbol{x}) - b = \sum_{i=1}^N c_i y_i k(\boldsymbol{x}_i, \boldsymbol{x}) - b \end{equation}

コラム:kernelの背後にある数学

カーネルで非線形変換後の内積を置き換えられる背景には,**再生核Hilbert空間(reproducing kernel Hilbert space)**の理論がある.本文では「変換後の座標を作らず,内積の値だけを直接求める」という計算上の働きをつかめれば十分だ.理論を詳しく学ぶ場合は,正定値kernel,Mercerの定理,再生核Hilbert空間を調べるとよい.


  1. 主問題 primal problem,双対問題 dual problem,強双対性 strong duality あたりをキーワードにして検索して勉強すると良い. ↩︎

レポート No.14

提出ファイルについての注意

近年はセキュリティ上の懸念から,実行形式のプログラムなどをメールに添付すると,受信側サーバがメールそのものを拒絶することがある.そういう問題を避けるため,レポートをファイルで提出するときは,実行形式などの危険視されやすい形式を避けよう.

要するに,レポートはPDFファイルにして送るのが良い,と思っておこう.

以下の課題について,自らの将来のスキルアップにつながるように調査と考察を行い,
     学籍番号-氏名-14.pdf
というファイルとしてレポートを作成し、 webフォーム から教官宛に提出しよう.

なお,レポートを $\TeX$ 等で作成したものを印刷した「紙媒体」を教官に直接手渡す形で提出してもよいが、物質によるレポート提出は常に破損や紛失の可能性があるのであまりお勧めはしないぞ.

課題

  1. SVM を適用すると良さそうな分類問題の例が LIBSVM Data に挙げられているので,この中からどれでもよいので対象として SVM を適用してみよう. やり方は今回の授業の最初の実習例のように自分で地道にプログラムしてから Optim package などで最適化をしても良いし,LIBSVM package を使っても良い.

  なお,LIBSVM Data のデータの取扱いに困るという人もいるだろう. たとえば,「データの意味がよくわからない」とか,「そもそもどうやって julia に取り込んだら良いのか技術的に困る」という理由でだ. そこでそのあたりを少し解説しておこう.

まず,データの意味について解説しておこう.といっても全部は多すぎるので,a1a から a9a という名前のものについてのみにしておこう. 他もこれで慣れればきっとだいじょうぶだろう.

まず,これら a1a - a9a のデータは,UCI にある Adult というデータを変換したもので,これはアメリカ人の特徴データ 14個と,それに加えて「年収が 5万ドル以上かどうか」を分類したデータだ.

たとえば a1a の 1行目は

-1 3:1 11:1 14:1 19:1 39:1 42:1 55:1 64:1 67:1 73:1 75:1 76:1 80:1 83:1

となっているが,最初の項 -1 が分類を表していて,この場合は年収が5万ドル以下ということになる. 5万ドルを越えるとここが +1 になるのだ.

そして残りの 14個のデータは「下記の特徴 14個を全部並べ直して 123個の特徴8に直し,そのうち,あてはまるものを列挙したもの」になる.

よって,例えば a1a の 1行目の第2項は 3:1 は年齢が下から 3番めの段階にあたること,第3項の 11:1 は雇用型が 6番目の型にあたること,を意味するのだ.

  • 14個の特徴とその並べ直しのための情報
仮番号 特徴 数字か種類か.数字の場合,何段階に変換したか.
01 年齢 数字 → 5段階に変換
02 雇用型 8種類
03 体重 数字 → 5段階に変換
04 学歴 16種類
05 教育暦年 数字 → 5段階に変換
06 結婚歴 7種類
07 職業 14種類
08 家族関係 6種類
09 人種 5種類
10 性別 2種類
11 資本利益 数字 → 2段階に変換
12 資本損失 数字 → 2段階に変換
13 週あたり労働時間 数字 → 5段階
14 出身国 41種類

注意: 上の合計をとって 5段階x3個 + 2段階x2個 + (8 + 16 + 7 + 14 + 6 + 5 + 2 + 41)種類 = 123 で,これをすべて「種類」とみなして 123種類になる.

そこでデータは 123次元の空間に属する点とみなして,該当する種類の次元の座標を 1 にして,そうでない次元の座標を 0 としているのだ. 当然,123次元のうち 14次元の座標が 1 で,残りの次元の座標は 0 になる.

なんでこんなことにしているかというと,データに数字のものと,属性タイプみたいなものが混在しているので,全体を一度に取り扱うためにはこうした方が良い,というもののようだ.

次にこのデータの読み込み方だが,少し注意すべきことが有る. データを先に少し眺めてみると気づくだろうが,欠損しているデータが混じっているのだ.

欠損データが混じっていると,読み込みだけでも面倒だし,読み込んだ後でも取扱いには注意が必要だ9. そこで,以下に示す「データの取扱いのためのDataFramespackage」の出番だ.

  データを扱うときは,欠損などがあっても問題なくデータを扱えるDataFrames package を使おう!

あとは小さな工夫をするだけなので,a1a データファイルを対象にした例として以下にプログラムを見せよう.

まず,データそのままだと扱いにくいので,CSV形式に直してファイルを作る.

 1# データを簡単に整形して CSV ファイルにする.
 2ifn = open("a1a", "r")
 3ofn = open("a1a.csv", "w")
 4
 5while !eof(ifn)
 6  line = readline(ifn) # 1行読み込む
 7
 8  line2 = replace( line,  s":1" => s"")
 9  line3 = replace( line2, r" $" => s"")
10  line4 = replace( line3, s" " => s",")
11    
12  println( ofn, line4 )
13end
14
15close(ifn)
16close(ofn)

次に,必要なパッケージを2つ読み込む. いずれも VIORAS 環境には既にインストール済みだ. 個人環境を使っている等の状況で未インストールならば先にインストールしておこう(インストール方法はもう省略する.いつもの通りだ).

1using CSV
2using DataFrames

次に,CSVファイルを読み込んで,データに取り込む.

1# CSV ファイルを読み込む.
2raw_data = CSV.File("a1a.csv", header = false)
3
4# データとして認識させる.
5data = DataFrame( raw_data )
1,605 rows × 15 columns (omitted printing of 6 columns)
	Column1	Column2	Column3	Column4	Column5	Column6	Column7	Column8	Column9
	Int64	Int64	Int64	Int64	Int64	Int64	Int64	Int64	Int64
1	-1	3	11	14	19	39	42	55	64
2	-1	3	6	17	27	35	40	57	63
…略…

次に,DataFrames package の機能をありがたく使わせてもらって,欠損データを排除しよう.

1# 欠損値があるデータを除去.これができるのが DataFrames のありがたさ.
2clean_data = dropmissing(data)

あとは読み込んだデータを SVM 用にするのだ.

 1data_num = size(clean_data)[1]
 2# きれいなデータの個数
 3
 4X_train = Matrix{Int}(undef, data_num, 123)
 5# 特徴データ用行列.属性 123個について,0 or 1 が入る.1 は14箇所に入る.
 6
 7Y_train = Vector{Int}(undef, data_num)
 8# 分類データ用ベクトル.-1 は収入が 5万ドル以下,+1 は5万ドルより大きいという意味.
 9
10for i in 1:data_num
11    
12    x = zeros(Int, 123) # 特徴データ用ベクトル.最初は全部 0 を入れておいて…
13    
14    y = ( clean_data[[i], 1] )[1] # 分類.DataFrames のデータの読み方は少しややこしいかもね.
15
16    for dn in 2:15 # 14個の特徴について
17        on_bit = ( clean_data[[i],dn] )[1]  # 何番目の属性が ON なのかを読んで,
18        x[on_bit] = 1 # その属性を +1 に設定する.
19    end
20    
21    X_train[i, :] .= x # 特徴データ用 Matrix に格納.
22    Y_train[i] = y     # 分類データ用 Vector に格納.
23end

こうして得られる X_train が学習用特徴データに,Y_train が学習用の分類データになる.中身を見てみると次のような感じになっているはずだ.

1X_train
2# これが学習用特徴データになる.
1478×123 Matrix{Int64}:
 0  0  1  0  0  0  0  0  0  0  1  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  1  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  0  1  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
…略…
1Y_train
2# これが学習用分類データになる.
1478-element Vector{Int64}:
 -1
 -1
 -1
 -1
…略…

  1. 二種類への分類は用途が限られすぎ…と思うかもしれないが,現在は多種類の分類機としての拡張もなされている.もちろん,二種類への分類を単純に繰り返しても多種類分類ができるぞ. ↩︎

  2. このことから,support vector machine の名前がある. ↩︎

  3. ソフトマージンを採用しないもの, たとえばハードマージンを採用する(サンプルデータ点をすべて分類成功するような超平面を求めることに相当する)ものももちろん SVM なのだが,よほどきれいなデータ点でないと「全体がそもそも失敗する」ので実用性が低くなってしまうな. ↩︎

  4. ベクトル $\boldsymbol{x}$ の二次形式とは,ある行列 $A$ が存在して $\boldsymbol{x}^{T} A \boldsymbol{x}$ と書けるものを言うぞ.要するにたくさんの変数 $\{x_i\}_{i=1}^n$ の二次式(ただし全部二次項)だな. ↩︎

  5. もちろん,「おおよそ状況がわかっている」ケースなどでは非線形変換のみでも問題ない. ↩︎

  6. カーネルを決めることが非線形変換を決めることにも相当するので,非線形変換を先に決める必要がないという言い方をしても良い. ↩︎

  7. Chih-Chung Chang and Chih-Jen Lin, LIBSVM : a library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2:27:1--27:27, 2011. Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm ↩︎

  8. 調べてみるとこれを 129個の特徴に変換して扱っている人もいる.そういう人は特徴の No.11,12(資本による利益値,損失値)をそれぞれ 5段階に変換しているようだ. ↩︎

  9. 欠損データをどうするかはもちろん方針次第だ.ただ,どうするにせよ,素でプログラムを書くと面倒な処理になりかねない. ↩︎