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$ はそれぞれこの超平面の法線ベクトルと切片で，超平面自身が

\[ \boldsymbol{w}\cdot\boldsymbol{x} - b = 0 \]

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

そして上の図で左上側にあるグループ「青」を $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}

と書き直せる．

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

\[ \frac{2}{\| \boldsymbol{w} \|} \]

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

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

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

すると式($\ref{eq:in-group}$) は \[ y_i (\boldsymbol{w}\cdot\boldsymbol{x}_i - b) \geq 1 \mbox{ for all } (\boldsymbol{x}_i, y_i) \label{eq:cond-hardmargin} \] と書き直せるので，ハードマージン下での問題全体は

  制約条件付き最小化問題:
      $\| \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) \]

とするのである． 条件を守れていない場合はこの値がプラスになることに留意しておこう．ちなみに条件を守れている場合はマイナスの値になる．

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

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

とすれば扱いやすいのでそうしよう．

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

\[ 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\} \]

を考えてこれを最小化すればいい，ということになる． つまり，

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

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

ちなみに，係数 $\lambda$ はペナルティをどれだけ軽視するかという量になっていて，大きくすると条件の守られ方が緩くなる，という関係にある． 逆に係数 $\lambda$ を小さくすると条件の守られ方は厳しくなり，ハードマージンに近くなる．

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

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

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

次に，ペナルティの計算と，超平面のパラメータの定数自由度 $C$ の計算の両方のための情報となる量 $y (\boldsymbol{w}\cdot\boldsymbol{x} - b)$ を計算する関数を作る． これは，条件を満たしていれば自由度 $C$ の手がかりに，満たしていなければペナルティに使われる情報なので，ここでわけておこう．

 1function r_check(x, y, w, b)
 2    r_max = 10000.0
 3    
 4  r =  y * ( w ⋅ x - b )
 5  # 条件を満たしていれば，r は正になる．
 6    
 7  if r > 0.0
 8        return r, 0.0
 9        # r_min候補値と，ペナルティを返す
10    else
11        return r_max, -r
12        # r_min候補値と，ペナルティを返す
13    end
14end

そして，データ全体に渡って計算して，$C$ と今の 自由定数$C$ を正しく調整してない $\boldsymbol{w}$, $b$ でのペナルティ値を計算する関数を用意する．

 1function condition_and_penalty(w,b)
 2  min_r = 10000.0
 3  ap = 0.0
 4  for i in 1:div(N,2)
 5
 6    res = r_check(Xp[:,i], 1.0, w, b)
 7    if res[1] < min_r 
 8      min_r = res[1]
 9    end
10    ap += res[2]
11
12    res = r_check(Xn[:,i], -1.0, w, b)
13    if res[1] < min_r 
14      min_r = res[1]
15    end
16    ap += res[2]
17
18  end
19
20  return 1/min_r, ap / N
21  # 1/min_r は C に, ap/N はペナルティ項に一致する．
22
23end

それから，自由定数 $C$ を用いて正しく最小化したい関数 $L$ を次のように定める．

1function func_L(v)
2    λ = 0.01
3    w = v[1:2]
4    b = v[3]
5    
6    C, penalty = condition_and_penalty(w, b)
7    
8    return λ*(C^2)*(norm(w)^2) + C*penalty
9end

あとは Optim package におまかせでいい．次のような感じだ． ちなみに探索初期値は，超平面が右上がりになっている，というぐらいの設定だ．

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

 * Status: success

 * Candidate solution
    Final objective value:     5.424325e-02

 * Found with
    Algorithm:     Nelder-Mead

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

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

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

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

([-11302.434571580347, 19945.169790257743], 26716.384970865736)

ただしこれは定数自由度 $C$ の調整を受けてない数字なので，正しく調整しておこう．

1C = condition_and_penalty( w, b )[1]
2w = C * w
3b = C * b

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

1w

2-element Vector{Float64}:
 -1.1482483276835298
  2.0262897972985368

1b

2.7141979164101984

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

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$ によって

\[ \boldsymbol{w} = \sum_{i=1}^N c_i y_i \boldsymbol{x}_i \]

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

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

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

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

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

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

非線形変換

まず非線形変換であるが，これは形式的には \[ \boldsymbol{x} \mapsto \boldsymbol{\varphi}(\boldsymbol{x}) \] とするだけだ． ただし，このあとカーネルトリックを採用するため，
この$\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}(\boldsymbol{x}_i)\cdot\boldsymbol{\varphi}(\boldsymbol{x}_j) \]

の項)を，カーネル(Kernel)と呼ばれる「計算量が少ない二変数関数」

\[ k( \boldsymbol{x}_i, \boldsymbol{x}_j ) \]

で置き換えてしまう操作を言うのだ． なにそれ？ と思うかもしれないが，再生核ヒルベルト空間 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 \\ & = & \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_{-}$ の元であると判定する．

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

カーネルの種類

ちなみにカーネルの種類等であるが，まあ典型的なのは次のような感じか．

• カーネルの例

  名前	関数 $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 を採用して実習してみよう． これをインストールしていない人はまずは以下のようにインストールしよう．

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

レポート

以下の課題について能う限り賢明な調査と考察を行い，
2023-AppliedMath7-Report-13
という題名をつけて e-mail にて教官宛にレポートとして提出せよ. なお，レポートを e-mail の代わりに TeX で作成した書面にて提出してもよい.

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

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

課題

1. 自分で適当な問題を作って，SVM を適用してみよう．
   
   
2. (チャレンジ問題) 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 を使おう！

あとは小さな工夫をするだけなので，順にプログラムを見せよう．
まず，データそのままだと扱いにくいので，以下のようにして 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)

次に，必要なパッケージの読み込みだ． なお，未インストールならば先にインストールしておこう．

1using CSV
2using DataFrames

次に，下記のようにして CSVファイルを読み込んで，データに取り込む． 欠陥データなどがあるので大量に Warning がでるが，気にしなくて良いだろう．

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,605個のデータが 1,478個に減っていて，約 130近くのデータは要素が揃ってないような「欠陥データ」だったことが見て取れる．

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

1,478 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
…略…

あとは読み込んだデータを 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
…略…