はじめに
こんにちは，(株)日立製作所 研究開発グループ サービスコンピューティング研究部の露木です。
機械学習の勉強を始めると「損失関数」や「正則化」といった言葉をよく見かけるかと思います。この正則化の効果を理屈と感覚で理解したいならば，シンプルな回帰モデルであるLasso回帰が題材として最適です。そこで本記事ではLasso回帰の特徴を簡単に説明した上で，scikit-learnのプログラムを動かしながらその振る舞いを調べます。特に，Lasso回帰のハイパーパラメータを変えながら非線形データの回帰を繰り返し，正則化の効果を確認します。
Lasso回帰の特徴
Lasso回帰は数値データに対して、線形回帰を行う手法の一つです。例えば特徴量が $n$ 次元あるときに，与えられたデータに対して以下の式の重み係数 $w_i$ と切片 $b$を最適化すること（線形回帰）を考えます。
$$ \hat{y} = \sum_{i=1}^n w_i x_i + b $$
例えば，この線形回帰のパラメータを最小二乗法で求める場合，MSE（Mean Squared Error, 平均二乗誤差) が最小となるような重み係数${\bf w}$と切片$b$を，適当な数 (例えば$m$個) の学習データから計算することになります。ここでMSEは
$$ {\rm MSE} = \frac{1}{m} \sum^m_{p=1} \left(y_p - \hat{y}_p \right)^2 $$
で計算できます。最小二乗法におけるMSEのような「学習する際に最小化する関数」を，機械学習の分野では損失関数(Loss function)と呼びます。
一方，Lasso回帰では正則化項として$w_i$のL1ノルム
$$ || {\bf w} ||_1 \equiv \sum_{i=1}^n | w_i | $$
をMSEに加えた，下記の値を損失関数に用います。ただし，$ \alpha $ はLasso回帰のハイパーパラメータであり，正則化の強さを制御します。
$$ {\rm MSE} + \alpha ||{\bf w}||_1$$
このような正則化項は重み係数の絶対値を小さくし，過学習を抑える効果があります。特に，L1ノルムはすべての$i$について $|w_i| = 0$ となるときに最小値をとるので，最小二乗法と比較してLasso回帰は係数（の一部）を0にし，予測に用いる特徴量の次元を削減できる特徴があります。
そのため，最小二乗法による線形回帰や，正則化項にL2ノルムを使うRidge回帰と比較して，Lasso回帰は分析データの特徴量の次元が高く，次元を削減したい（不要な特徴量が多く含まれる可能性が高い）場合に使用されることが多い手法です。
実行例と解説
以下では，実際にLasso回帰を使ってサンプルの数値データを学習し，回帰モデルを作成します。ソースコードはPython3のJupyterノートブックで実行する前提であることに注意してください。
ライブラリの準備
まず必要なライブラリがインストールされていない場合は，次のセルをコメントアウトしてから実行してインストールしてください。
# !pip install numpy scikit-learn matplotlib pandas
次にライブラリを読み込み，初期設定を行います。
# ライブラリを読み込む
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Lasso, LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
# Jupyterでグラフ表示するためのおまじない
%matplotlib inline
データ準備
本記事ではサンプルデータとして，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("x1")
plt.ylabel("y")
plt.scatter(X, y, color="black")
生成したサンプルデータの70%はモデルの学習に使う学習データとし，30%はモデルの精度評価に使うテストデータとします。モデルは，そもそも学習データを再現するように学習するものですから，新しいデータに対するモデルの精度（汎化性能）を正しく評価するには学習データとは別のテストデータを使わなければなりません。
# ソートが楽になるようにdataframe形式にする
X = pd.DataFrame(X)
y = pd.DataFrame(y)
# 学習データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=2018)
# データがランダムソートされているので直す
X_train = X_train.sort_index()
y_train = y_train.sort_index()
X_test = X_test.sort_index()
y_test = y_test.sort_index()
# 学習データとテストデータのグラフを表示
plt.scatter(X_test, y_test, label="test data")
plt.scatter(X_train, y_train, label="training data")
plt.title("sample data")
plt.xlabel("x1")
plt.ylabel("y")
plt.legend()
Lasso回帰の学習
次にLasso回帰を学習します。Lasso回帰のハイパーパラメータのうち，重要なのは$\alpha$です。最初に述べたように，Lasso回帰では下記の式を最小化するモデルを学習します。
$$ \frac{1}{m} \sum^m_{p=1} \left(y_p - \hat{y}_p \right)^2 + \alpha ||{\bf w}||_1$$
上記の式を見ればわかるように，$\alpha$は正則化項$||{\bf w}||_1$の係数であり，$\alpha$が大きいほど正則化の効果が強くなります。一般に最適な$\alpha$の値は，$\alpha$の値を変えながら学習と精度評価を繰り返しながら精度最大となる値を探すこと（ハイパーパラメータ探索）で決定しますが，今回はとりあえず $\alpha = 0.05$ とし、他のパラメータはデフォルトのまま学習します。
なお，その他のパラメータについてはscikit-learnの公式ドキュメントを参照下さい。
clf = Lasso(alpha=0.05)
clf.fit(X_train, y_train)
学習したモデルの切片$b$は次のように確認できます。
print("切片b = %f" %(clf.intercept_))
切片b = 0.168758
次に，学習したモデルの重み係数${\bf w}$を確認します。今回の場合，生成したサンプルデータの特徴量は1次元なので重み係数も一次元になります。
print("重み係数w")
print(clf.coef_)
重み係数w
[0.18507533]
モデルを使った予測
今度は学習したモデルを使ってテストデータの性質を正しく再現できるか確認します。以下のグラフを見ると，Lasso回帰は線形回帰ですからテストデータの右肩上がりな性質は再現できていますが，非線形的なsinc関数の振る舞いは再現できないことがわかります。
# 予測値の計算
p = clf.predict(X_test)
# グラフ化
plt.scatter(X_test, y_test,label="test data")
plt.title("predicted results")
plt.xlabel("x1")
plt.ylabel("y")
plt.plot(X_test, p, label="model")
plt.legend()
精度の検証
最後に回帰モデルの精度として平均二乗誤差の値と，決定係数 $R^2$ の値を確認します。決定係数は1に近いほど精度良くデータを再現していることを示します。先程のグラフでも明らかですが， $R^2$ 値からも今回学習したモデルの精度はとても低いと言えます。
print("平均二乗誤差 MSE = %f" % (mean_squared_error(y_test, p)))
print("決定係数 R^2 = %f" % (r2_score(y_test, p)))
平均二乗誤差 MSE = 0.197404
決定係数 R^2 = 0.389922
カーネル法による非線形データへの対応
回帰モデルの精度が低い理由は，今回のサンプルデータが非線形データであるにも関わらず，強引に線形回帰をしているからです。このような場合は，まずは適当なカーネル関数 $k({\bf x}_p, {\bf x}_q)$ を用いた非線形的な写像で特徴量空間の基底を張りなおしてから，改めて線形回帰して非線形データへ対応します。ここで，$p$と$q$はデータ点を示すラベルです。
このテクニックはカーネル法やカーネルトリックなどと呼ばれ，Support Vector Machineをはじめとした機械学習の分野で広く使われています。カーネル関数の具体的な形式としてはガウス関数が頻繁に用いられます。すなわち，
k({\bf x}_p, {\bf x}_q) = \exp(-\gamma ||{\bf x}_p - {\bf x}_q||^2)
です。なお，$\gamma$はハイパーパラメータとなる係数です。このようにガウス関数を利用したカーネル関数は，2つの入力ベクトル ${\bf x}_p $ と ${\bf x}_q$ 間の距離に依存するため，RBFカーネル (Radial Basis Function カーネル，動径基底関数) とも呼ばれます。カーネル法の理論的な背景はカーネル多変量解析やパターン認識と機械学習 (下) といった教科書を参照いただくとして，ここでは計算の手続きのみ紹介します。
例えば，簡単のために2次元 $(n=2)$ の特徴量空間において ${\bf x}_1 = (1, 2)$ となるデータ点1 と ${\bf x}_2 = (3, 4)$となるデータ点2 からなる学習データ ${\bf X} = ({\bf x}_q, {\bf x}_p)$ があったとします。この${\bf X}$ にカーネル法を適用する場合，カーネル関数を行列要素とするグラム行列 ${\bf K}$を計算します。
\begin{eqnarray}
{\bf K} = \left(
\begin{array}{ccc}
k({\bf x}_1, {\bf x}_1) & k({\bf x}_1, {\bf x}_2) \\
k({\bf x}_2, {\bf x}_1) & k({\bf x}_2, {\bf x}_2)
\end{array}
\right)
\end{eqnarray}
Scikit-Learnのデフォルト値は $\gamma = 1/n$ であり，今考えている学習データの数は $n=2$ ですからRBFカーネルの値は
k({\bf x}_p, {\bf x}_p) = \exp (- \frac{1}{2} \times ||(1, 2) - (1, 2)||^2 ) = 1
k({\bf x}_p, {\bf x}_q) = \exp (- \frac{1}{2} \times ||(1, 2) - (3, 4)||^2 ) = 0.018
k({\bf x}_q, {\bf x}_q) = \exp (- \frac{1}{2} \times ||(3, 4) - (3, 4)||^2 ) = 1
となります。つまり，グラム行列は
\begin{eqnarray}
{\bf K} = \left(
\begin{array}{ccc}
1 & 0.018 \\
0.018 & 1
\end{array}
\right)
\end{eqnarray}
です。カーネル法では，このような計算から作成したグラム行列の1行目 $(1, 0.018)$ と2行目 $ (0.018, 1)$を，${\bf x}_1$や${\bf x}_2$に代わる新たな特徴量として線形回帰を学習することになります。
実際に，先程のサンプルデータでカーネル法とLasso回帰の組み合わせを試してみると，以下のグラフのように精度良く非線形データを再現できることがわかります。
# RBFカーネルを用いて、データを非線形写像する
kX_train = rbf_kernel(X_train, X_train)
# 写像したデータを用いて、Lasso回帰を学習する
clf = Lasso(alpha=0.01, max_iter=2000)
clf.fit(kX_train, y_train)
# グラフ描画用途に，細かいグリッド点で計算
kX_plot = rbf_kernel(X_plot, X_train)
p_plot = clf.predict(kX_plot)
# グラフ化
plt.scatter(X_test, y_test, label="test data")
plt.title("predicted results")
plt.xlabel("x1")
plt.ylabel("y")
plt.plot(X_plot, p_plot, label="model")
plt.legend()
精度の検証
次にカーネル法を使って学習したモデルの精度を計算します。学習データと同じ方法でテストデータの特徴量もカーネル法を適用しなければならないことに注意してください。例えば，$ {\bf x}_3 $ がテストデータであった場合，新たな特徴量 $(k({\bf x}_3, {\bf x}_1), k({\bf x}_3, {\bf x}_2)) $ を計算してから，モデルに入力します。
実際に，平均二乗誤差の値と決定係数 $R^2$ の値を確認すると，カーネル法によって決定係数は1に近づいてモデルの精度が改善したことがわかります。
# 予測値の計算
kX_test = rbf_kernel(X_test, X_train)
p = clf.predict(kX_test)
print("平均二乗誤差 MSE = %f" % (mean_squared_error(y_test, p)))
print("決定係数 R^2 = %f" % (r2_score(y_test, p)))
平均二乗誤差 MSE = 0.075589
決定係数 R^2 = 0.766393
正則化項の効果
Lasso回帰の切片$b$と重み係数${\bf w}$を確認すると，カーネル法によって次のように変化しています。もともと生成したサンプルデータの特徴量は1次元でしたが，グラム行列の次元，つまりは学習データの数 $n$ まで重み係数の次元が増加しています。ただし，多くの場合は $w_i=0$となっており，Lasso回帰の正則化項によって適切に次元数が抑えられたことがわかります。
print("切片: %f" %(clf.intercept_))
print("重み係数")
print(clf.coef_)
切片: -0.320667
重み係数
[-0. -0. -0. -0. -0.02787288 -0.
-0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0.21635477
-0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0.
-0. -0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0.1345954 0.85714895
0.25103181 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0.35264187 0.53731577 0.06222342 0.
0. 0. ]
正則化をしない場合
もし正則化をしない場合は，どのような問題が発生するのでしょうか。これを確認するために，カーネル法と最小二乗法で線形回帰を実施してみます。
以下の結果を見れば分かるように，正則化しなかった結果として10の12乗を超えるような非常に大きい重み係数 $w_i$ が学習されてしまいます。
# カーネル法で写像したデータを用いて、最小二乗法で回帰を学習する
clf = LinearRegression()
clf.fit(kX_train, y_train)
print("切片: %f" %(clf.intercept_))
print("重み係数w")
print(clf.coef_)
切片: 23309336.993793
重み係数w
[[-1.06175813e+11 -2.77640256e+12 1.26765873e+13 -5.12623261e+12
-1.38273791e+13 8.79604003e+11 7.69837966e+12 4.74701548e+12
2.11218344e+12 -7.46896027e+11 8.85306741e+10 -6.83877865e+12
-4.34070599e+12 -3.76833474e+12 1.55829107e+12 2.57514527e+12
4.33176053e+12 3.42120342e+12 2.43700236e+12 1.75191078e+12
5.53895415e+11 -5.73364576e+12 -3.48700973e+12 -4.61993566e+12
1.51489024e+12 7.16450476e+12 8.53118228e+11 6.86393427e+12
-2.61838471e+12 -3.15174892e+12 -3.43594518e+12 -4.09724668e+12
-5.08240220e+12 -1.48092049e+12 1.49466382e+12 3.54513253e+12
7.78206474e+12 1.35350194e+12 5.28751146e+11 -2.94718437e+12
-4.49940866e+12 -3.26032749e+12 -7.52061817e+11 2.74156489e+12
3.74941589e+12 2.00718096e+12 2.34428690e+12 1.90774598e+12
-1.77347088e+11 -2.02046805e+12 -6.92936084e+12 -2.81387214e+12
-2.49413931e+12 3.28036180e+12 2.17334102e+12 5.65521153e+12
3.02428364e+12 1.40609675e+11 -1.12198198e+12 -1.81170834e+12
-3.21986026e+12 -1.70136773e+12 -3.79928373e+12 -5.86705813e+11
-7.53575429e+11 1.63883856e+12 6.22817182e+11 1.66293297e+12
2.93554120e+12 3.04040172e+12 3.43134874e+12 -3.93912018e+11
-8.55312025e+11 -1.38545774e+12 -2.95846285e+12 -3.79637943e+12
6.22907801e+11 -2.85748679e+12 -7.66826946e+11 3.28779085e+12
3.33565504e+12 2.56240457e+12 3.41778160e+12 1.70530134e+12
1.41319483e+12 -2.31418039e+12 -3.44909140e+12 -3.49971728e+12
-3.31989914e+12 1.08443009e+11 -2.23746259e+12 -7.60936917e+11
1.08551996e+12 -1.09931627e+12 2.78151024e+12 1.00133548e+12
2.72349583e+12 1.94583200e+12 2.72402529e+12 2.06330753e+12
-2.27935291e+12 -3.45253950e+11 -3.68320093e+12 -2.79696805e+12
-2.55690293e+12 -2.31223689e+12 2.70056909e+12 1.59286185e+12
2.56090675e+12 3.91329213e+12 9.53151868e+11 1.23552771e+12
-1.57644338e+12 -2.30290399e+12 -1.18371258e+11 -3.43042704e+12
-2.50579329e+11 -3.21307682e+12 -1.13671895e+12 1.84580356e+12
2.23469233e+12 1.38535332e+12 2.93220673e+12 1.35744664e+12
-1.50729149e+11 1.34962673e+12 -3.13414414e+12 -4.36742444e+12
-1.34411190e+11 -4.87516738e+11 1.35300867e+12 2.81225427e+12
1.20272696e+12 1.07918522e+12 -1.85124722e+12 -2.76806707e+12
-1.76479225e+12 1.91973212e+12 1.88091974e+12 -8.93844040e+11]]
このモデルによる予測値を確認すると，次のグラフのように過学習 (overfitting)していることがわかります。正則化をしないとパラメータ数が増えすぎて，複雑すぎるモデルを学習してしまうといえます。
# plot用の予測値を計算
kX_plot = rbf_kernel(X_plot, X_train)
p_plot = clf.predict(kX_plot)
# グラフ化
plt.scatter(X_test, y_test, label="test data")
plt.title("predicted results")
plt.xlabel("x1")
plt.ylabel("y")
plt.ylim(-2, 2)
plt.plot(X_plot, p_plot, label="model")
plt.legend()
また，モデルの精度を確認しても，正則化しない場合は極端に悪化すると言えます。これは，前述のグラフ右端の点において非常に大きな値をとっていることに由来します。
# 予測値の計算
kX_test = rbf_kernel(X_test, X_train)
p_test = clf.predict(kX_test)
print("平均二乗誤差 MSE = %f" % (mean_squared_error(y_test, p_test)))
print("決定係数 R^2 = %f" % (r2_score(y_test, p_test)))
平均二乗誤差 MSE = 4.380572
決定係数 R^2 = -12.538149
最後に
本稿ではLasso回帰の特徴を説明した上で，実際に正則化の役割を数値実験で示しました。線形モデルを非線形データへ対応可能とする代わり特徴量の次元が増えるカーネル法を利用する場合，Lasso回帰の正則化によって不要な次元を削減することが過学習回避に効果的といえます。
なお，本稿ではわかりやすさのためにLasso回帰を題材としましたが，なにかの問題に実適用する場合は より性能の良い カーネルリッジ回帰 や サポートベクトル回帰 の採用を検討してください。手前味噌ではありますが， サポートベクトル回帰の説明は過去に執筆しております ので，そちらも合わせてご覧いただければ幸いです。