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}

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

つまり, 条件を満たす範囲でなるべく $\| \boldsymbol{w} \|$ を小さくすればよいという，最小化問題なのだ．

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

すると式($\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} と書き直せるので，ハードマージン下での問題全体は

  制約条件付き最小化問題:
      $\| \boldsymbol{w} \|$ を最小化せよ．
      ただし，全てのデータ $(\boldsymbol{x}_i, y_i)$ に対して $y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 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)$ であることを考えると，

$1 - y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b)$ がマイナスならば条件を守れていて，プラスならば守れていない

ということになる． しかも，守れていない場合はこの量の大きさそのものが「どれくらい条件を守れていないのか」という「大きさ，程度」も表してくれそうだ．

そこで，この条件に関するペナルティ $P$ を

\begin{equation} P(\boldsymbol{x}_i, y_i; \boldsymbol{w}, b) := \left\{ \begin{array}{lcl} 1 - y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) & : & \mbox{左の量がプラスの場合} \cr 0 & : & \mbox{otherwise} \end{array}\right. \end{equation}

とすれば今まで考えていることと整合する． というわけで，ペナルティをこう設定しよう．

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

\begin{equation} L(\boldsymbol{w}, b; \lambda) := \lambda \| \boldsymbol{w} \|^2 + \frac{1}{N} \sum_{i=1}^N P(\boldsymbol{x}_i, y_i; \boldsymbol{w}, b ) \label{eq:softmargin-loss} \end{equation}

に対し，$\boldsymbol{w}, b$ を調整して最小化すればいい，ということになる．

  上の式 (\ref{eq:softmargin-loss}) を最小化すると $y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 1$ かもしれないが サポートベクター上で$y_i(\boldsymbol{w}\cdot\boldsymbol{x}_i - b) > 1$ となっていたりしないか？ と思うかもしれない． 実はその心配は要らないのだ． というのも，式 (\ref{eq:softmargin-loss}) の最小化では第一項の $\| \boldsymbol{w} \|^2$ を小さくしていくときに第二項でペナルティをギリギリまで小さくできるところまでいくわけだが，その限界がまさにサポートベクター上での条件 $y_i(\boldsymbol{w}\cdot\boldsymbol{x}_i - b) = 1$ を満たすときだからだ． まあ，もっときちんと数式で納得したい人は，サポートベクター上で $y_i(\boldsymbol{w}\cdot\boldsymbol{x}_i - b) = 1$ となるという条件無しでの議論を一回真面目にやってみると良い．

さて話を戻すと，つまりは

  最小化問題:
      与えられた正係数 $\lambda$, データ $(\boldsymbol{x}_i, y_i)_i^N$ に対し， 量 $L$ を最小化する $\boldsymbol{w}, b$ を求めよ．

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

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

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

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

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

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

個人環境で Optim が未インストールの場合は下記のようにしてインストールしておこう．

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,:])

データ(点配置と分類)をまとめておこう．

1Xall = hcat( Xp, Xn )
2Labels = vcat( [ 1.0 for i in 1:div(N, 2) ],  [ -1.0 for i in div(N,2)+1: N ] )

次に，ペナルティを計算する関数を作る．

 1function all_penalty(w, b)
 2    p = 0.0
 3    n = length(Labels)
 4
 5    for i in 1:n
 6      r = 1.0 - Labels[i] * ( w ⋅ Xall[:, i] - b )
 7      p += (r > 0.0 ? r  : 0.0)
 8    end
 9
10    return p / n
11end

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

それから，最小化したい肝心の関数 $L$ を次のように定める．

1function func_L(v, λ)
2    w = v[1:2]
3    b = v[3]
4
5    return λ*(norm(w)^2) + all_penalty( w, b )
6end

あとは Optim package におまかせでいい．次のような感じだ． 計算はほぼ一瞬で終わる．

ちなみに, $\lambda = 0.0001$ は何度か試しながら適当に設定した値で，探索初期値は「超平面が右上がりになっている」というぐらいの適当な設定だ．

1res = optimize( v -> func_L(v, 0.0001), [-1.0, 1.0, 1.0 ])
2# 小さくしたい関数と，探索初期値を渡す．

 * Status: success

 * Candidate solution
    Final objective value:     5.381904e-04

 * Found with
    Algorithm:     Nelder-Mead

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

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

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

1v = res.minimizer
2w, b = v[1:2], v[3]

([-0.8369112692899874, 2.167152048840457], 4.31619962132352)

つまり，超平面のパラメータは次のような値に近い値になるはずだ(データを乱数で作っているので，やる人ごとに少し違う)．

1w

2-element Vector{Float64}:
 -0.8369112692899874
  2.167152048840457

1b

4.31619962132352

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

1scatter(Xp[1,:], Xp[2,:])
2scatter!(Xn[1,:], Xn[2,:])
3
4plot!( x -> ( -w[1]*x + b)/w[2], label = "hyperplane" )

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

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

1classifier(x) = sign( w⋅x - b)

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

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$ についてはハードマージンのときと同じだ．

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

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

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

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

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

非線形変換

まず非線形変換であるが，これは形式的には \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 パッケージのインストール

個人環境で未インストールの場合は下記のようにしてインストールしておこう．

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"],
Int32[1, 2, 3], Float64[], Int32[],
LIBSVM.SupportVectors{Vector{InlineStrings.String15}, Matrix{Float64}}(30,
…略…
LIBSVM.SVMNode(1, 6.0), LIBSVM.SVMNode(1, 5.8), LIBSVM.SVMNode(1, 6.3)]),
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 は良い性能を持つのだな，と実感できる．

レポート

下記要領でレポートを出してみよう．

• e-mail にて，
• 題名を 2024-numerical-analysis-report-13 として,
• 教官宛(アドレスは web の "TOP" を見よう)に，
• 自分の学籍番号と名前を必ず書き込んで，
• 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，

下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．

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

まあ要するに，レポートは 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
…略…