正則化の意味を理解する

「過学習」に陥らないようにモデルを構築することは、機械学習（教師有り学習）の基本です。正則化は過学習を防ぐテクニックの一つです。Weight Decayとも呼ばれ、DeepLearningなど最新の手法でも使われています。

正則化は難しい概念ではないので、実際に使ってみるとわかります。下記のサンプルコードを動作させるには、

Scikit-learn
NumPy
Matplotlib

の3つのライブラリが必要です。

今回は線形回帰を例に、Scikit-learnのライブラリを使って、正則化の挙動を確認します。その前にデータを用意しましょう。今回は描画を行うために、２次元のデータセットで確認します。

但し、一般には特徴量が多い（描画できない）場合に、過学習が起きやすい傾向があることに留意しましょう。

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot') 

rng = np.random.RandomState(1)
x = 10 * rng.rand(50)
y = np.sin(x) + 0.1 * rng.randn(50) 

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x, y)
ax.set_xlabel('x')
ax.set_ylabel('y')

plt.show()

f:id:roughnessjoker0815:20190331025232p:plain

このトイデータセットを使ってxを説明変数、yを目的変数とする回帰を行います。ここで、過学習が起こる条件の1つは、特徴量が多くモデルの表現力が高いことです。

今回は、xについて基底関数を適用することでモデルを複雑さ（表現力）を意図的に増やして、過学習を起こします。例えば、xに関する高次の項を用意するために、scikit-learnのPolynomialFeatures変換器が利用できます。

$y= a_0 + a_1x + a_2x^ 2 + a_3x^ 3$

のように、xの3次の項まで用意したい場合は、PolynomialFeaturesのdegreeオプションに3を指定してインスタンスを生成すれば良いです。

poly = PolynomialFeatures(degree=3)
# 切片項が不要な場合、include_bias=Falseを指定する。
# poly = PolynomialFeatures(degree=3, include_bias=False)  
x_deg_3 = poly.fit_transform(x[:, np.newaxis])

今回はガウス基底関数の和をモデルにすることで、より顕著な過学習を確認することが出来るため、ガウス基底関数を利用します。ガウス基底関数はscikit-learnに組み込まれていないので、自前で用意する必要があります。（ガウス分布＝正規分布のことです。xの区間をN等分して、N個の中心を持つ正規分布の和をモデルとします。）

下記のコードは、書籍「Pythonデータサイエンスハンドブック」の引用です。

from sklearn.base import BaseEstimator, TransformerMixin

class GaussianFeatures(BaseEstimator, TransformerMixin):
    #1次元入力の等間隔ガウス関数
    
    def __init__(self, N, width_factor=2.0):
        self.N = N
        self.width_factor = width_factor
        
    @staticmethod
    def _gauss_basis(x, y, width, axis=None):
        arg = (x - y) / width
        return np.exp(-0.5 * np.sum(arg ** 2, axis))
    
    def fit(self, X, y=None):
        # データ範囲に沿ってN個の中心を作成する
        self.centers_ = np.linspace(X.min(), X.max(), self.N)
        self.width_ = self.width_factor * (self.centers_[1] - self.centers_[0])
        return self
    
    def transform(self, X):
        return self._gauss_basis(X[:, :, np.newaxis], self.centers_,
                                 self.width_, axis=1)

正則化を行わずに、30次のガウス基底関数の和を使って、回帰を行ってみましょう。

from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LinearRegression

model = make_pipeline(GaussianFeatures(30),
                      LinearRegression())
model.fit(x[:, np.newaxis], y)
xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用
y_pred = model.predict(xfit[:, np.newaxis])

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x, y)
ax.plot(xfit, y_pred)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax2 = fig.add_subplot(212)
ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_)
ax2.set_xlabel('basis')
ax2.set_ylabel('coefficient')

plt.show()

f:id:roughnessjoker0815:20190331032154p:plain

上が回帰プロット、下が30次のガウス基底それぞれの係数プロットです。確かに全ての散布図データにフィットするような回帰曲線が引かれていますが、未知のデータに対応できるモデルとは言えないでしょう。機械学習（特に教師有り学習）では、汎化性能のみに興味があり、既に得られているデータへの当てはまりは全く重要でない場合が多いというところがポイントです。これが、過学習が起きている状態です。

係数プロットからはかなり大きい係数の値を取っていることが読み取れます。

このように複雑なモデルでは、大きな係数を取ることで学習データに対してのみ上手くフィットする極端な形のグラフがかけてしまうのです。

そこで、係数にペナルティをかけて、極端な形のグラフをかけないようにしようというのが正則化のアイデアです。

数式的には、損失関数に正規化項を加えます。

正則化なし: 損失関数 $L = \sum_{i=1}^n (y_i - \hat y_i)^ 2$
正則化あり: 損失関数 $L = \sum_{i=1}^n (y_i - \hat y_i)^ 2 +$ （正則化項）

正則化項には、L2正則化であれば、

$\lambda \sum_{j=1}^ m w_j^ 2$

L1正則化であれば、

$\lambda \sum_{j=1}^ m |w_j|$

が使われます。

回帰係数は、損失関数を最小にするように求めるため、残差を小さくする他に、正則化項も小さくすること、つまり回帰係数を小さくすることが要求されるのです。

主要な線形回帰モデル

scikit-learnの主要な線形回帰器は4つあります。正則化を行わないLinearRegression、L2正則化を行うRidge、L1正則化を行うLasso、L2正則化とL1正則化のハイブリッドElasticNetです。それぞれ表にまとめます。

モデル	正則化	正則化パラメータ（デフォルト）
LinearRegression	正則化なし	なし（正則化ありモデルのalpha=0に相当）
Ridge	L2正則化	alpha=1
Lasso	L1正則化	alpha=1
ElasticNet	L1 + L2正則化	alpha=1, l1_ratio=0.5

Ridge

Ridge回帰（L2正則化）は、最も一般的な正則化の形式です。回帰係数の二乗和（L2ノルム）を正則化項に加えることで、係数を抑えます。

この記事では行いませんが、全ての正則化に共通することとして、正則化の強さを決めるハイパーパラメータを最適化することが重要です。

from sklearn.linear_model import Ridge

model = make_pipeline(GaussianFeatures(30),
                      Ridge(alpha=0.1))
model.fit(x[:, np.newaxis], y)
xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用
y_pred = model.predict(xfit[:, np.newaxis])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(211)
ax.scatter(x, y)
ax.plot(xfit, y_pred)
ax.scatter(6, model.predict(np.array([[6]])), marker='*', s=200)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax2 = fig.add_subplot(212)
ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_)
ax2.set_xlabel('basis')
ax2.set_ylabel('coefficient')

plt.show()

f:id:roughnessjoker0815:20190331041848p:plain

正則化項を追加することで、係数のスケールが大幅に小さくなり、感覚的に正しいグラフが得られました。

Lasso

Lasso回帰（L1正則化）は、回帰係数の絶対値の和（L1ノルム）を正則化項に加えることで、係数を抑えます。

L1正則化の重要な特性として、不要な特徴量の回帰係数を（正則化の強さ次第で）0にできることが挙げられます。そのため、特徴選択の用途で使われることがあります。

from sklearn.linear_model import Lasso

model = make_pipeline(GaussianFeatures(30),
                      Lasso(alpha=0.001, max_iter=100000))
model.fit(x[:, np.newaxis], y)
xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用
y_pred = model.predict(xfit[:, np.newaxis])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(211)
ax.scatter(x, y)
ax.plot(xfit, y_pred)
ax.scatter(6, model.predict(np.array([[6]])), marker='*', s=200)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax2 = fig.add_subplot(212)
ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_)
ax2.set_xlabel('basis')
ax2.set_ylabel('coefficient')

plt.show()

f:id:roughnessjoker0815:20190331043140p:plain

いくつかの係数が、完全にゼロになっていることが確認できます。

ElasticNet

Elastic Netは、L1正則化とL2正則化の両方を使用することで、より良い正則化を目指したモデルです。

Elastic Netには正則化の強さを決めるパラメータが2つあり、1つは全体としての正則化の強さを決めるalpha、もう1つはL1正則化とL2正則化の強さの比率を決めるl1_ratioがあります。

ハイパーパラメータを上手く調整すれば、Ridge回帰やLasso回帰よりも高い性能が得られることがある一方、適切な正則化パラメータを見つける労力がかかるというデメリットもあります。

from sklearn.linear_model import ElasticNet

model = make_pipeline(GaussianFeatures(30),
                      ElasticNet(alpha=0.001, l1_ratio=0.5, max_iter=100000))
model.fit(x[:, np.newaxis], y)
xfit = np.linspace(0, 10, 1000) # 回帰直線プロット用
y_pred = model.predict(xfit[:, np.newaxis])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(211)
ax.scatter(x, y)
ax.plot(xfit, y_pred)
ax.scatter(6, model.predict(np.array([[6]])), marker='*', s=200)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax2 = fig.add_subplot(212)
ax2.plot(np.linspace(0, 10, 30), model.steps[1][1].coef_)
ax2.set_xlabel('basis')
ax2.set_ylabel('coefficient')

plt.show()

f:id:roughnessjoker0815:20190331044316p:plain

Lasso回帰と比較して、いくつかの特徴量の係数が0ではなくなっていることがわかります。

まとめ

今回は復習も兼ねて、正則化についてまとめました。実際のデータ分析で過学習が起きる状況は、特徴量が多数あり、少なくともそのままでは回帰の様子を可視化できないことがほとんどだと思います。そのため、２次で過学習のイメージを持っておくことは、正則化の直感的な理解に役立つと思います。

この記事では回帰の例を使って説明しましたが、ロジスティック回帰などの線形分類の正則化項についても同様の考え方で理解することができます。ただし、ロジスティック回帰の正則化パラメータのCは、 $\lambda$ の逆数であるので、Cを大きくするほど正則化が弱まることには注意が必要です。

（線形サポートベクターマシンについては、Cはソフトマージンのパラメータなので、正則化とは異なりますが、モデルのバイアス・バリアンスを制御するという意味では同じです。）

ディープラーニングについても、損失関数を最小化するように係数（重み）を更新するという点については同じであり、全ての重みからなるノルムを損失関数に加えることで、重みにペナルティをかけて過学習を防ぐことを期待しています。

では、おやすみなさい🐼