Photo by Craig Whitehead on Unsplash

Suppor Vector Machine (SVM) とは

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

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

1. 空間上に既知のデータ点を配置する. このとき，データ点がおおまかに2つに別れて配置されている，と想定する．
   
   
2. 一つの超平面を用意し「その両側で空間を二分して」，そのどちらにデータ点があるかでデータを 2つの種類にうまく分類できるようにする．
   言い換えると，そのような適切な超平面を既知のデータ点から計算する．これが SVM の学習過程である．
     このとき，この超平面に一番近いデータ点を Support Vector ²と呼ぶ.
   
   
3. 未知のデータ点に対し，上で求めた超平面のどちらにあるかで，分類を判断する．

重要なポイントをおおまかに先に書いておく

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

1. 疑念.     空間にデータ点を配置しても，超平面でうまく二分できるようにデータ点が配置されているとは限らないのでは？
   
   回答.     そのとおり．そのため，通常はデータ点を
   
      ・非線形変換する
      ・(上の変換後のデータと組み合わせて)データ点の配置空間の次元を上げる
   
   ことでうまく二分できる「可能性を上げる」． これについては， SVM with polynomial kernel visualization (HD) (Youtube 動画, by udiprod)が大変わかりやすいイメージ動画になっているので一度見ておこう．
   さらに，
   
      ・うまく分離できないデータ点の存在を許すソフトマージンを採用する³
   
   ことで全体の柔軟性を上げるのだ． より詳しいことは後述する．
   
   
2. 疑念.     「適切な」超平面を計算する手順は？ とくに上記の非線形変換などを行うと，うまくできる気がしないが…
   
   回答.     「なるべく良い」超平面を探す最適化問題に帰着させることとちょっとした工夫で，既存の知識で可能になる． 重要なポイントを書いておこう．
   
      ・分類失敗するデータ点も許す「ソフトマージン」で超平面を計算する！
   
      ・非線形変換を行う場合，最適化問題の主問題ではなく双対問題を解く
   
      ・この双対問題はデータ点の二次形式⁴になる
   
      ・次元やデータ点数が多くなると二次形式の計算量増大がつらい
   
      ・二次形式相当の結果をいきなり計算できるカーネルを採用することで，計算量増大をゆるやかにできる．これはカーネルトリックと呼ばれる重要な工夫だ．
   

詳しくは後述する.

データ点をそのまま配置してほぼ分離が可能な場合

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

ハードマージンを考える: 全てのサンプルデータをきれいに分離する超平面へ

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

ただし，$\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} \|$ にペナルティを足して小さくすべき全体量を定義する

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

より詳しく解説しよう．

まず，本来守るべき制約条件 $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$ を小さくするとペナルティが重視され，マージンをハードマージンにとった状況に近くなる．

  実際にやってみる: その 1

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

なお，最適化のために 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,:])

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

そのために，まずは \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 = の式の右辺で w と Xp などの間にある記号 ⋅ (見えにくいが)は, ベクトルの内積を意味する記号で，\cdot & tabキー で入力できる． dot 関数 の sytax 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" )

ふむ，たいへんうまく分離できているようだね． これで 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 による分類器が正しく作れた，ということになるね． 簡単な問題の場合はこんな感じだ．

さて，話の続きに戻ろう． ソフトマージンの考え方の続きだ．

双対問題

そもそも最適化問題の理論に，もとの問題(主問題 primary problem と呼ぶ)に対して，双対問題 dual problem と呼ばれる問題を考えると，

  双対定理(duality theorem):
      主問題か双対問題が最適解を持つならば，もう一方も最適解を持ち，かつそれらは一致する．

という定理がある． であるので，主問題と双対問題の「難易度」が異なる場合，易しい方を解けば良い，というアイディアにつながるわけだ．

詳細⁵は省くが 上のソフトマージンの最小化問題にも双対問題を考えることができ，しかもこちらの方が話を拡張しやすいのだ． そこでこれを見ておこう．以下のようになる．

  上の最小化問題の双対問題:
      与えられた正係数 $\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$ についてはハードマージンのときと同じだ．

データ点をそのまま配置すると分離が難しい場合(非線形変換へ)

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

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

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

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

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

非線形変換

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

双対問題

次に双対問題を考える点についてだが，これは下記のように 式$(\ref{eq:normal-dual-problem})$を少し変形した問題になる．

  非線形変換時の双対問題(カーネルトリック採用前):
      与えられた正係数 $\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}

の項)の形で現れることに着目して，この項を $\boldsymbol{\varphi}$ の内積ではなく， カーネル(Kernel)と呼ばれる「計算量が少ない二変数関数」

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

で置き換えてしまう操作を言うのだ． なにそれ？ と思うかもしれないが， 再生核ヒルベルト空間 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}$ は計算しなくてよい(やっぱり $\boldsymbol{\varphi}$ を必要としないのだ)． 解説しておくと，次のように変形を考えれば，カーネルだけあれば良いことがわかる．

• 超平面の切片 $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}

カーネルの種類

ちなみにカーネルの種類等であるが，まあ典型的なのは次のような感じだ． よくわかんないときはとりあえず 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$

  実際にやってみる: その 2

ここでは SVM を扱うちょっと有名なライブラリ LIBSVM⁸ を 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 ] # 検証用

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

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

で， 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 は良い性能を持つのだな，と実感できる．

レポート No.13

以下の課題について能う限り賢明な調査と考察を行い，
     学籍番号-氏名-13.pdf
というファイルとしてレポートを作成し、 webフォーム から教官宛に提出しよう.

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

  注意
　　近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする． そういうことを避けるため，レポートをファイルで提出するときはそういった懸念のあるファイル形式のものではないようにしよう．

まあ要するに，レポートは pdf ファイルにして送るのが良い ということだと思っておこう．

課題

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

  なお，LISVM 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個の特徴⁹に直し，そのうち，あてはまるものを列挙したもの」になる．

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

• 14個の特徴とその並べ直しのための情報

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

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

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

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

欠損データが混じっていると，読み込みだけでも面倒だし，読み込んだ後でも取扱いには注意が必要だ¹⁰． そこで，以下に示す「データの取扱いのための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
…略…