1. はじめに
こんにちは，(株)日立製作所 研究開発グループ サービスコンピューティング研究部の露木です。
分類問題に使うサポートベクトルマシン (SVM) は有名ですが，これを数値データの回帰予測に応用したアルゴリズムとして SVR (Support Vector Regression, サポートベクトル回帰) があります。今回は，SVRのハイパーパラメータの役割を理解した上で，設定できるようになることを目指します。
本記事の構成は，まず2節でSVRの定義に基づいてハイパーパラメータの役割の理論的な解釈を行います。次に，3節ではsinc関数にノイズを加えた人工データについて，実際にハイパーパラメータを変えながら学習を行い，ハイパーパラメータの影響を感覚的にも理解できるようにします。
なお，本記事ではハイパーパラメータを変えながら結果を確認する「機械学習の実験管理」を一つのノートブックの中で行いますが，このアプローチは実験数や機械学習チームの規模が大きくなると破綻します。明日の記事では，同じSVRのハイパーパラメータ選択を題材としながらも，オープンソースの実験管理ソフト MLflow を利用してシステム化する方法を紹介します。ぜひそちらも合わせてご覧ください。
2. 理論
最初に定義式に基づいてSVRのハイパーパラメータの役割を説明します。そのような説明が不要な場合は，次節の実行例からご覧ください。
SVMはもともと2値分類問題を解くために開発されたアルゴリズムであり，カーネル関数の柔軟性から分類問題では定番といえるほど広く使われているアルゴリズムです。分類問題にSVMを適用する場合は，特にSVC (Support Vector Classification，サポートベクトル分類)と呼ばれます。一方，回帰問題にも対応できるように，目的変数を連続値へ拡張したSVRもあります。SVRは非線形な回帰問題を比較的精度良く解けるため，SVCほどではありませんがSVRも広く使われています。
SVRを使うにあたって設定するハイパーパラメータは主に3種類あり，正則化係数 $C$ と不感度係数 $\varepsilon$ とカーネル関数です。各ハイパーパラメータの意味を理解するために，まずは線形回帰を行う単純なSVRモデル $f({\bf \it x}_i)$ を考えます。
$$
f({\bf \it x}_i) = {\bf \it w}^{\rm T} {\bf \it x}_i + b \tag{1}
$$
$\bf \it x_i$ は説明変数であり， $\bf \it w$と$b$は学習時に最適化する重みベクトルとバイアスです。$\bf \it w$と$b$を決めるSVRの最適化問題は下記の式(2)で表現され，ここにハイパーパラメータである $C$ と $\varepsilon$ が現れます。式(2)の第一項は正則化項であり，重み係数の絶対値を小さくすることで過剰適合 (Overfitting) を回避します。第二項は損失関数であり，目的変数の真値 $y_i$ と予測値 $({\bf w}^{\rm T} x_i +b )$の差，つまり残差 (予測誤差) の総和 $\sum _i | y_i - ({\bf w}^{\rm T} x_i +b )|$を最小化する$\bf \it w$と$\bf\it b$を決定するための項です。
$$
({\bf \it w},b) = \arg \min_{{\bf w},b} \frac{1}{2} ||{\bf \it w} ||^2 + C \sum_{i} \max [ | y_i - ({\bf w}^{\rm T} x_i +b ) | - \varepsilon, 0 ] \tag{2}
$$
式(2)の第二項に含まれる$\max [| y_i - ({\bf w}^{\rm T} x_i +b ) | - \varepsilon,0 ]$ は，残差の絶対値 $|y_i - ({\bf w}^{\rm T} x_i +b ) |$ が $\varepsilon$ 以下であれば0を用いることを意味します。つまり，$\varepsilon$の役割は訓練時に真値と予測値の誤差が小さいサンプルを無視することであり，このような小さなノイズを無視して頑健なモデルを学習することだといえます。SVMの用語では，学習に影響するサンプルを「サポートベクトル」と呼びますが，$\varepsilon$が大きいほど無視するサンプルの数が増加し，サポートベクトルの数は減少します。
一方，式(2)より，$C$は正則化項と損失関数のどちらを重視して最適化するか決めるハイパーパラメータであることがわかります。$C$が大きいほど，訓練データに対する残差を重視することになるため，大きすぎる$C$は過剰適合を招きます。
ここまでは式(1)のように線形回帰を行うSVRを説明しました。カーネル関数の選択は非線形の回帰問題へSVRを拡張するために，式(1)の内積を (マーサーの定理を満たす) 任意の関数で置き換える操作を意味します。カーネル関数の導入によって，特徴量ベクトル ${\bf \it x}_i$ の基底を非線係な変換で張り直し，新たに得た特徴量ベクトル${\bf \it \phi}({\bf \it x}_i)$についてモデルパラメータ${\bf \it w}$と$b$を最適化することで，最適化問題そのものは線形回帰と同一でありながらも非線形的な回帰問題へ対応可能になります。
カーネル関数の具体的な表式は多数ありますが，多くの場合で精度が良く，頻繁に用いられるのはガウス関数を使ったRBFカーネル (Radial Basis Function Kernel)です。その他，scikit-learnでは多項式近似とシグモイド関数を選択できます。適切なカーネル関数はデータの形状によって異なりますが，まずRBFカーネルを試し，精度が不十分であれば交差検証に基づいて精度の良い関数を選択するべきです。
本節の内容について詳細は，下記の参考文献をご覧ください。
- 竹内一郎，鳥山昌幸，機械学習プロフェッショナルシリーズ サポートベクトルマシン (2015)，講談社 https://www.kspub.co.jp/book/detail/1529069.html
次節の実行例では，実際に$C$と$\varepsilon$，そしてカーネル関数を変えながら予測結果を確認していきます。
3. 実行例
3.1. 環境準備
※ ここより先のコードはPython3であり，Jupyterノートブックでの実行を想定しています。
まずは，必要なライブラリをインストールします。今回はscikit-learnの実装を利用してSVRを学習します。
# ライブラリがインストールされていない場合のみ，このセルをコメントアウトして実行する
# !pip install numpy scikit-learn matplotlib
ライブラリをインポートします。
# ライブラリを読み込む
import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVR
from sklearn.model_selection import train_test_split
# Jupyterノートブックでグラフを表示するためのおまじない
%matplotlib inline
3.2. データ準備
訓練データとして，sinc関数に右肩上がりな成分とノイズを加えたデータを生成します。
# サンプルデータ数
m = 200
# 乱数のシード値を指定することで，再現性を保つ
np.random.seed(seed=2018)
# 「-3」から「3」の間で等間隔にm個のデータを作成
X = np.linspace(-3, 3, m)
# 後のグラフ描画用途に，10倍細かいグリッドを準備しておく
X_plot = np.linspace(-3, 3, m*10)
# 周期的なsin関数(第一項)に右上がり成分(第二項)と乱数(第三項)を加えたデータを作る
y = np.sinc(X) + 0.2 * X + 0.3 * np.random.randn(m)
# グラフ表示するため，各数列を1列の行列に変換
X = X.reshape(-1, 1)
y = y.reshape(-1, 1)
X_plot = X_plot.reshape(-1,1)
# グラフを表示
plt.title("sample data")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.scatter(X, y, color="black")
生成したサンプルデータの70%はモデルの学習に使う「訓練データ」とし，30%はモデルの精度評価に使う「テストデータ」とします。モデルは，そもそも訓練データを再現するように学習するものですから，未知データに対するモデルの精度（汎化性能）を正しく評価するには訓練データに含まれないデータ（テストデータ）を使わなければなりません。
# 訓練データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=2019)
# 訓練データとテストデータのグラフを表示
plt.scatter(X_test, y_test, label="test data", edgecolor='k',facecolor='w')
plt.scatter(X_train, y_train, label="training data", color='c')
plt.title("sample data")
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend()
3.3. SVRモデルの学習
まずは，3種類のカーネル関数 (linear, ploy, rbf) で学習を行います。linearは線形回帰，ployは多項式展開，rbfはガウス関数が使用されます。
reg_linear = SVR(kernel='linear', C=1, epsilon=0.1, gamma='auto')
reg_poly = SVR(kernel='poly', C=1, epsilon=0.1, gamma='auto')
reg_rbf = SVR(kernel='rbf', C=1, epsilon=0.1, gamma='auto')
reg_linear.fit(X_train, np.ravel(y_train))
reg_poly.fit(X_train, np.ravel(y_train))
reg_rbf.fit(X_train, np.ravel(y_train))
[セルの出力]
SVR(C=1, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',
kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False)
精度指標として決定係数$R^2$を比較します。$R^2$は1に近いほど，良くデータを近似できていることを示します。
scores = (reg_linear.score(X_test, y_test),
reg_poly.score(X_test, y_test),
reg_rbf.score(X_test, y_test))
plt.bar(("Linear", "poly", "RBF"), scores)
plt.xlabel("Kernel")
plt.ylabel("$R^2$ score")
plt.show()
ここではRBFカーネルが最も精度が良い結果となりました。また，多項式(poly)カーネルを用いると，線形(linear)カーネルよりも精度が悪くなっています。非線形データを対象にしているにもかかわらず，線形カーネルのほうが高精度であることは不自然です。次節で可視化しながら，その理由を調べます。
3.4. 予測結果の確認とハイパーパラメータの調整
決定係数$R^2$だけではどのようなモデルを得られたのか不明であるので，予測値を可視化して調べてみます。まず，予測値をグラフ表示するための関数を定義します。
# 描画用の関数を定義
def plot_result(model, X_train, y_train, score):
# 予測値の計算
p = model.predict(np.sort(X_test))
# グラフ化
# plt.scatter(X_test, y_test,label="test data")
plt.clf()
plt.scatter(X_test, y_test, label="test data", edgecolor='k',facecolor='w')
plt.scatter(X_train, y_train, label="Other training data", facecolor="r", marker='x')
plt.scatter(X_train[model.support_], y_train[model.support_], label="Support vectors", color='c')
plt.title("predicted results")
plt.xlabel("$x$")
plt.ylabel("$y$")
x = np.reshape(np.arange(-3,3,0.01), (-1, 1))
plt.plot(x, model.predict(x), label="model ($R^2=%1.3f$)" % (score), color='b')
plt.legend()
(a) カーネルによる予測結果の変化
次に線形カーネルの予測値を図示します。図中では青色の直線がモデルによる予測値を示し，データ点は訓練データとテストデータを表します。また，赤いx印は$\varepsilon$の範囲内であり，学習に影響しなかったデータを表します。
図を見ると線形カーネルはテストデータの右肩上がりな性質は再現できていますが，非線形的なsinc関数の振る舞いは再現できないことがわかります。当然ではありますが，線形回帰は非線形なデータに適していないといえます。
plot_result(reg_linear, X_train, y_train, reg_linear.score(X_test, y_test))
次に多項式カーネルの結果を確認します。多項式カーネルによってモデルが曲線を描けるようにはなっていますが，scikit-learnのデフォルト設定では3次多項式で展開するため，今回のデータに対しては表現力が不足していると考えられます。
plot_result(reg_poly, X_train, y_train, reg_poly.score(X_test, y_test))
そこで，多項式カーネルの次元数を6次元まで増やして再度学習します。degree=6 が次元数の設定です。
reg_poly = SVR(kernel='poly', C=1, epsilon=0.1, gamma='auto', degree=6, coef0=1)
reg_poly.fit(X_train, np.ravel(y_train))
[セルの出力]
SVR(C=1, cache_size=200, coef0=1, degree=6, epsilon=0.1, gamma='auto',
kernel='poly', max_iter=-1, shrinking=True, tol=0.001, verbose=False)
多項式カーネル (degree=6) の精度を他の場合と比較すると，RBFカーネルと同等程度まで向上しました。
scores = (reg_linear.score(X_test, y_test),
reg_poly.score(X_test, y_test),
reg_rbf.score(X_test, y_test))
plt.bar(("Linear", "poly (degree=6)", "RBF"), scores)
plt.xlabel("Kernel")
plt.ylabel("$R^2$ score")
plt.show()
6次元多項式のカーネルによる予測結果を可視化します。やはり，データをよく近似できる曲線が得られていることがわかります。
plot_result(reg_poly, X_train, y_train, reg_poly.score(X_test, y_test))
最後にRBFカーネルの予測結果を確認します。多項式カーネルとの差はわずかですが，両端において不自然な下降がなく，より正しく近似できていることがわかります。一般にRBFカーネルは計算時間が早く精度も良い場合が多いため，特別な理由が無い限りはRBFカーネルを選ぶべきです。
plot_result(reg_rbf, X_train, y_train, reg_poly.score(X_test, y_test))
(b) 不感度係数 ε の効果
次に，不感度係数$\varepsilon$の効果を確かめます。RBFカーネルで，$\varepsilon$の値を0，0.1，0.9とした場合の3つのモデルを学習します。
model1 = SVR(kernel='rbf', C=1, epsilon=0.0,gamma='auto').fit(X_train, np.ravel(y_train))
model2 = SVR(kernel='rbf', C=1, epsilon=0.1, gamma='auto').fit(X_train, np.ravel(y_train))
model3 = SVR(kernel='rbf', C=1, epsilon=0.9, gamma='auto').fit(X_train, np.ravel(y_train))
3つのモデルの精度 ($R^2$値) を比較します。$\varepsilon =0$に対して $\varepsilon =0.1$ のモデルはわずかに改善していますが， $\varepsilon =0.9$ のモデルは大幅に精度が悪化していることがわかります。
scores = (model1.score(X_test, y_test), model2.score(X_test, y_test), model3.score(X_test, y_test))
plt.bar(("0.0", "0.1", "0.9"), scores)
plt.xlabel("$\\varepsilon$")
plt.ylabel("$R^2$ score")
plt.show()
$\varepsilon = 0$の場合は，すべての訓練データが学習に影響するサポートベクトルとなります。
plot_result(model1, X_train, y_train, model1.score(X_test, y_test))
$\varepsilon = 0.1$の場合は，モデル近傍のデータ点 (赤いx印)がサポートベクトルから取り除かれます。誤差の小さいノイズを無視することで，テストデータに対する精度 ($R^2$値) がわずかに改善していることがわかります。
plot_result(model2, X_train, y_train, model2.score(X_test, y_test))
一方，$\varepsilon = 0.9$の場合は，値が大きすぎるために訓練データの大半 (赤いx印) がサポートベクトルから取り除かれてしまっています。サポートベクトルの数が少なすぎる結果として，テストデータに対する精度 ($R^2$値) が大幅に悪化した事がわかります。
plot_result(model3, X_train, y_train, model3.score(X_test, y_test))
以上のように $\varepsilon$ はゼロでは最適と言えませんが，大きすぎると精度を極めて悪化させるといえます。最適な$\varepsilon$ の値はデータに応じて変化しますが，ばらつきに比例して最適な$\varepsilon$も大きくなります。
(c) 正則化係数Cの効果
最後に，正則化係数$C$の効果を確かめます。RBFカーネルで，$C$の値を0.01，1，1000とした場合の3つのモデルを学習します。
model1 = SVR(kernel='rbf', C=0.01, epsilon=0.1, gamma='auto').fit(X_train, np.ravel(y_train))
model2 = SVR(kernel='rbf', C=1, epsilon=0.1, gamma='auto').fit(X_train, np.ravel(y_train))
model3 = SVR(kernel='rbf', C=1000, epsilon=0.1, gamma='auto').fit(X_train, np.ravel(y_train))
3つのモデルの精度 ($R^2$値)を比較すると$C=1$の場合がバランスが良く，高精度になっていることがわかります。
scores = (model1.score(X_test, y_test), model2.score(X_test, y_test), model3.score(X_test, y_test))
plt.bar(("0.01", "1", "1000"), scores)
plt.xlabel("$C$")
plt.ylabel("$R^2$ score")
plt.show()
モデルの予測値を確認すると，$C=0.01$の場合は平坦な予測曲線となり，テストデータを全く再現できていないことがわかります。この場合，正則化項を過剰に重視した結果，予測誤差を小さくするパラメータを最適化できていないといえます。
plot_result(model1, X_train, y_train, model1.score(X_test, y_test))
$C=1$の場合は，正則化項と損失関数のバランスが良く，過剰適合せずに予測誤差を小さくするパラメータを最適化できたことがわかります。
plot_result(model2, X_train, y_train, model2.score(X_test, y_test))
一方，$C=1000$の場合は損失関数を過剰に重視し，訓練データに対して過剰適合していることがわかります。
plot_result(model3, X_train, y_train, model3.score(X_test, y_test))
以上のように $C$ も，小さすぎても大きすぎても精度を悪化させるハイパーパラメータです。最適な$C$ の値はデータに応じて変化しますので，精度などの指標を確認しながら選択することが大切です。
最後に
本記事ではSVRのハイパーパラメータを理解した上で設定できるようになることを目指し，理論と数値実験の両面から説明を行いました。SVR利用時に設定すべきハイパーパラメータは3種類あり，正則化係数 $C$ と不感度係数 $\varepsilon$ とカーネル関数です。
まず，$C$は過剰適合を防ぐ正則化の強さを制御するパラメータです。 $C$ が小さいほど正則化が強くなり，平坦な予測曲線を描くようになります。次に，$\varepsilon$ は小さな誤差を含むサンプルを無視し，ノイズに対して頑健なモデルを学習するための不感度係数です。 $\varepsilon$ が大きすぎると最適化に必要なサンプルまで無視してしまい，正しく学習できなくなってしまいます。最後に，カーネル関数は，SVRを非線形最適化問題へ拡張する関数です。まずは，一般に広く用いられるRBFカーネルを試し，精度が不十分であれば交差検証に基づいて精度の良い関数を選択するべきといえます。
以上の内容に関して，本記事では1つのノートブックでハイパーパラメータを変えながら数値実験を繰り返すアプローチをとりました。このようなアプローチでは，実験数や開発メンバーが増えると破綻します。明日の記事では，同じSVRのハイパーパラメータ選択を題材としながらも，オープンソースの実験管理ソフト MLflow を利用してシステム的に管理する方法を紹介します。ぜひそちらも合わせてご覧ください。