C-SVMのチューニング(正則化パラメータの最適化)

# はじめに
C-SVMのチューニングでは，正則化パラメータ($C$)の最適値をCross Validationで求める方法が提案されています．
本エントリーでは，pythonのscikit-learnを使って，SVMのクラス分類に必要な$C$を最適化する方法を紹介します.
SVMの損失関数は，データ数に依存して変化するため，データ数のスケーリングによって，$C$の最適値を評価する必要があります．

# 正則化(Regularization)
SVMで分類を行う際，データが線形分離できない場合，分離超平面のマージン内のデータを許す，ソフトマージンSVMを用います．
その際，分類しないデータについてペナルティを与え，ペナルティの重要度を正則化パラメータ($C$)によって決定します．
$C→∞$とすると，誤分類を許容しないハードマージンSVMとなります．
ソフトマージンSVMでは，以下の式で表される誤差を最小とする最適化を行います．

```math
C \sum_{i=1}^n L \left (  f ( x_i ), y_i  \right ) + \Omega ( \omega ) 
```
　$C$:正則化パラメータ
　$L$:損失関数
　$Ω$:ペナルティ関数(正則項)
　$n$:データ数
　$f$:判別関数
　$x_i$:入力変数
　$y_i$:クラスラベル{$+1,-1$}

正則化によって，過学習を抑えて，汎化能力を高めることができます.
上記の式では，$C$，$Ω$を追加項として与えることで，正則化を行います．
$C$は$L$と$Ω$の寄与率を制御するパラメータとなります．
$L$，$Ω$には，以下の関数のいずれかを選択します．

　$L$:hinge関数，あるいはhinge関数の二乗(squared_hinge)

```math
hinge(m) = max \{ 1-m, 0 \}
```
```math
squared\_hinge(m) = ( max \{ 1-m, 0 \} )^2
```

　$Ω$:$L1$正則化($p=1$)，あるいは$L2$正則化($p=2$)

```math
\Omega ( \omega ) = \dfrac {1}{p} \Vert f \Vert^p
```

# scikit-learnを用いたC-SVM
pythonの機械学習ライブラリ[scikit-learn](http://scikit-learn.org/stable/index.html)では,線形分離でSVMのクラス分類を行う[LinearSVC](http://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVC.html#sklearn.svm.LinearSVC.score)があります.
LinearSVCでは,$L$,$Ω$をそれぞれ設定することができます.

> Similar to SVC with parameter kernel='linear', but implemented in terms of liblinear rather than libsvm, so it has more flexibility in the choice of penalties and loss functions and should scale better to large numbers of samples.
> https://github.com/scikit-learn/scikit-learn/blob/51a765a/sklearn/svm/classes.py#L18-L21

一方，scikit-learnにある[SVC](http://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html#sklearn.svm.SVC)では,$L=$hinge,$Ω=L2$で変更はできません.

> loss : string, 'hinge' or 'squared_hinge' (default='squared_hinge')
> Specifies the loss function. 'hinge' is the standard SVM loss (used e.g. by the SVC class) while 'squared_hinge' is the square of the hinge loss.
> penalty : string, 'l1' or 'l2' (default='l2')
> Specifies the norm used in the penalization. The 'l2' penalty is the standard used in SVC. The 'l1' leads to ``coef_`` vectors that are sparse.
> https://github.com/scikit-learn/scikit-learn/blob/51a765a/sklearn/svm/classes.py#L33-L41

# C-SVMのチューニング
$C$の最適値をCross Validationによって求めることで,C-SVMをチューニングします．
但し，$L$はデータ数に応じて増加するため,$C$をデータ数に応じてスケーリングします．

> If we consider the loss function to be the individual error per sample, then the data-fit term, or the sum of the error for each sample, will increase as we add more samples. The penalization term, however, will not increase.
> http://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVC.html#sklearn.svm.LinearSVC

## ソースコード
ソースコードはscikit-learnにある"Scaling the regularization parameter for SVCs"の[Example](http://scikit-learn.org/stable/auto_examples/svm/plot_svm_scale_c.html)を元に作成しました.

```py:plot_svm_scale_c.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import sklearn.svm as svm

from sklearn.svm import LinearSVC
from sklearn.cross_validation import ShuffleSplit
from sklearn.grid_search import GridSearchCV
from sklearn import datasets

# データセットからirisを入手
iris = datasets.load_iris()
X, Y = iris.data, iris.target
# 4個の入力変数のうち、最初の2変数のみ使用
X = X[:,0:2]
n_samples = Y.size

# 正則化項のペナルティ関数、損失関数
# (penalty, loss) = (l2, l1), (l1, l2), (l2, l2)
# パラメータCの値域をnp.logspaceで設定
clf_sets = [(svm.LinearSVC(penalty='l2', loss='hinge', dual=True, tol=1e-3), np.logspace(-1, 2, 30), X, Y),
            (svm.LinearSVC(penalty='l1', loss='squared_hinge', dual=False, tol=1e-3), np.logspace(-1, 2, 30), X, Y),
            (svm.LinearSVC(penalty='l2', loss='squared_hinge', dual=True, tol=1e-3), np.logspace(-1, 2, 30), X, Y)]

for fignum, (clf, cs, X, Y) in enumerate(clf_sets):
    plt.figure(fignum, figsize=(9, 10))

    print 'penalty=%s, loss=%s' % (clf.penalty, clf.loss)
    print 'fraction, C, C*n_samples, max_score'
    for k,train_size in enumerate(np.linspace(0.3, 0.7, 3)[::-1]):
        param_grid = dict(C=cs)
        grid = GridSearchCV(clf, refit=False, param_grid=param_grid,
                            cv=ShuffleSplit(n=n_samples, train_size=train_size, n_iter=250, random_state=1))
        grid.fit(X, Y)
        scores = [x[1] for x in grid.grid_scores_]
        # cross validationで使用したデータ数,最大正答率でのパラメータC,C(スケーリング後),最大正答率
        print '%4.2f, %8.3f, %8.3f, %4.3f'%(train_size, cs[np.argmax(scores)], cs[np.argmax(scores)] * n_samples * train_size, np.max(scores))
        # パラメータCをサンプル数でスケーリングしない場合、スケーリングする場合のグラフ作成
        scales = [(1, 'No scaling'),((n_samples * train_size), '1/n_samples'),]

        for subplotnum, (scaler, name) in enumerate(scales):
            plt.subplot(2, 1, subplotnum + 1)
            plt.xlabel('C')
            plt.ylabel('CV Score')
            grid_cs = cs * float(scaler)  # scale the C's
            plt.semilogx(grid_cs, scores, label="fraction %3.2f" % (train_size))
            plt.title('scaling=%s, penalty=%s, loss=%s' % (name, clf.penalty, clf.loss))

    plt.legend(loc="best")
plt.show()
```

## 出力結果
- $Ω=L2$,$L=hinge$

| 訓練データ数(比率) | $C$ | $C*$訓練データ数 | 最大正答率 |
|:----------:|:----:|:------------:|:---------:|
| 0.70       |   11.721 | 1230.707 |  0.758 |
| 0.50       |   11.721 |  879.077 | 0.737 |
| 0.30       |   14.874 | 669.308 | 0.720 |

![iris_P2_L1_dual.png](https://qiita-image-store.s3.amazonaws.com/0/63906/c82a287b-ec66-0cb3-80f5-3ae3d2062aa0.png)

- $Ω=L1$,$L=squared\\_hinge$

| 訓練データ数(比率) | $C$ | $C*$訓練データ数 | 最大正答率 |
|:----------:|:----:|:------------:|:---------:|
|0.70 |   62.102| 6520.678| 0.779|
|0.50 |   62.102| 4657.627| 0.775|
|0.30 |   23.950| 1077.762| 0.767|

![iris_P1_L2_prime.png](https://qiita-image-store.s3.amazonaws.com/0/63906/05a61b2f-aaec-8ad3-097e-2faa0a20184f.png)

- $Ω=L2$,$L=squared\\_hinge$

| 訓練データ数(比率) | $C$ | $C*$訓練データ数 | 最大正答率 |
|:----------:|:----:|:------------:|:---------:|
|0.70 |   2.807|  294.758| 0.772 |
|0.50 |   5.736|  430.211| 0.766 |
|0.30 |   4.520|  203.416| 0.751 |

![iris_P2_L2_dual.png](https://qiita-image-store.s3.amazonaws.com/0/63906/d1182f30-ed18-51e4-5533-10a32e863d37.png)

## 考察

理論的には,$Ω=L1$では,データ数に依存するbiasがあるため,$C$をデータ数でスケーリングすることで性能を評価します.
一方,$Ω=L2$は,データ数に依存しないため,$C$のスケーリングなしでも評価が可能です.

> **l1-penalty case**
> In the l1 case, theory says that prediction consistency (i.e. that under given hypothesis, the estimator learned predicts as well as a model knowing the true distribution) is not possible because of the bias of the l1. It does say, however, that model consistency, in terms of finding the right set of non-zero parameters as well as their signs, can be achieved by scaling C1.

> **l2-penalty case**
> The theory says that in order to achieve prediction consistency, the penalty parameter should be kept constant as the number of samples grow.

> http://scikit-learn.org/stable/auto_examples/svm/plot_svm_scale_c.html

理論通り，$Ω=L1$,$L=squared\\_hinge$の結果は，データ数でスケーリングした$C$を用いることで，データ数に依存しないプロットが得られました．
また，$Ω=L2$,$L=squared\\_hinge$，$Ω=L2$,$L=squared\\_hinge$の結果についても，$C$のスケーリングしなくとも，データ数に依存しないプロットが得られ，理論と一致しました．
いずれの結果も，スケーリングすることで結果が明瞭になるので，C-SVMのチューニングには，データ数でスケーリングした$C$を用いたグラフで実施すると良いと言えます．

# チューニング前後のC-SVM結果
最後に，$C$をデフォルト値(1.0)，最適値の平均値を用いた，C-SVMの結果を示します．

## ソースコード
ソースコードは，先ほどの*plot_svm_scale_c.py*の*for fignum...*以下を次の様に修正しました.

```py:plot_separating_hyperplane.py
for fignum, (clf, cs, X, Y) in enumerate(clf_sets):
    plt.figure(fignum,figsize=(9,10))
    
    max_cs = []
    for k,train_size in enumerate(np.linspace(0.3, 0.7, 3)[::-1]):
        param_grid = dict(C=cs)
        grid = GridSearchCV(clf, refit=False, param_grid=param_grid,
                            cv=ShuffleSplit(n=n_samples, train_size=train_size, n_iter=250, random_state=1))
        grid.fit(X, Y)
        scores = [x[1] for x in grid.grid_scores_]

	max_cs.append(cs[np.argmax(scores)] * n_samples * train_size)

    x_min = 4.0
    x_max = 8.0
    y_min = 1.5
    y_max = 4.5

    XX, YY = np.mgrid[x_min-0.5:x_max+0.5:200j, y_min-0.5:y_max+0.5:200j]

    print 'penalty=%s, loss=%s' % (clf.penalty, clf.loss)
    print 'c, predict(average)'

    sets = [(1.0, 'default'),(np.average(max_cs) / Y.size, 'optimization'),]
    predicts = []
    for subplotnum, (c, name) in enumerate(sets):
	clf.set_params(C=c)
	clf.fit(X, Y)
	Z = clf.predict(np.c_[XX.ravel(), YY.ravel()])
	Z = Z.reshape(XX.shape)
	plt.subplot(2, 1, subplotnum + 1)
	plt.xlabel('Sepal length')
	plt.ylabel('Sepal width') 
	plt.xlim(x_min, x_max)
	plt.ylim(y_min, y_max)
 	plt.scatter(X[:, 0], X[:, 1], c=Y, zorder=10, cmap=plt.cm.Paired)
 	plt.axis('tight')
	plt.pcolormesh(XX, YY, Z, cmap=plt.cm.Paired)
        plt.title('C=%.3f, penalty=%s, loss=%s' % (c, clf.penalty, clf.loss))
	print '%.3f, %6.3f' % (c, clf.score(X,Y))

plt.show()
```

## 出力結果
- $Ω=L1$,$L=squared\\_hinge$

| $C$ | 正答率(平均値) |
|:----------:|:----:|
|1.000 |  0.720 |
|6.176|  0.787 |

![svm_p_l2_l_l1.png](https://qiita-image-store.s3.amazonaws.com/0/63906/4b188966-b69a-0870-55c6-f23c95f387a6.png)

- $Ω=L1$,$L=squared\\_hinge$

| $C$ | 正答率(平均値) |
|:----------:|:----:|
| 1.000 |  0.807 |
| 21.463 |  0.807 |

![svm_p_l1_l_l2.png](https://qiita-image-store.s3.amazonaws.com/0/63906/37085798-4592-e0e6-fea7-37dbb22f69d2.png)

- $Ω=L2$,$L=squared\\_hinge$

| $C$ | 正答率(平均値) |
|:----------:|:----:|
| 1.000 |  0.787 |
| 1.982 |  0.807 |

![svm_p_l2_l_l2.png](https://qiita-image-store.s3.amazonaws.com/0/63906/17fe125e-f882-e334-e083-36eb67e64f1d.png)

# 参考
[金森敬文「統計的学習理論」(講談社、2015年)](http://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E7%9A%84%E5%AD%A6%E7%BF%92%E7%90%86%E8%AB%96-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%83%97%E3%83%AD%E3%83%95%E3%82%A7%E3%83%83%E3%82%B7%E3%83%A7%E3%83%8A%E3%83%AB%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA-%E9%87%91%E6%A3%AE-%E6%95%AC%E6%96%87/dp/4061529056)
[Scaling the regularization parameter for SVCs　(scikit-learn公式ページ)](http://scikit-learn.org/stable/auto_examples/svm/plot_svm_scale_c.html)