正則化最小二乗法 - 機械学習に詳しくなりたいブログ

訓練データ数に対して多項式の次元数が大きすぎると、過学習が発生することを以前確認しました。（参考：線形回帰を最小二乗法で解く）過学習が発生するとき、係数 $\mathbf w$ が大きな値をとる傾向があるようです。よって、係数を小さい値に制限することができれば過学習が抑止できるということになります。この考えに基づいた手法が正則化です。最小二乗法の解の導出の式(5)で示した二乗和誤差

$E(\mathbf{w}) = \displaystyle \frac{1}{2} \|\mathbf X \mathbf w - \mathbf t \|^2 \tag{1}$

に対して、係数が大きくなることに対してペナルティを与える正則化項（罰金項）を以下のように追加します。

$E(\mathbf{w}) = \displaystyle \frac{1}{2} \|\mathbf X \mathbf w - \mathbf t \|^2 + \frac{\lambda}{2}\| \mathbf{w}\|^2 \tag{2}$

$\mathbf w$ が大きくなれば、そのぶん誤差関数も大きい値をとってしまうということですね。ここで $\lambda$ は正則化係数で、解析者が任意に決定する値です。式(2)を最小化する $\mathbf w$ を求めるには、最小二乗法の解の導出で導出したのと同様に展開して $\mathbf w$ で微分します。すると以下の式が求められます。

$\mathbf{w} = (\mathbf{X}^{T}\mathbf{X} + \lambda \mathbf{I})^{-1} \mathbf{X}^{T}\mathbf{t} \tag{3}$

導出の過程はほとんど同じなので省略します。僕は $\mathbf{X}^{T}\mathbf{X}\mathbf{w} + \lambda\mathbf{w} = (\mathbf{X}^{T}\mathbf{X} + \lambda \mathbf{I})\mathbf{w}$ の変形がちょっと戸惑いました。（ $\mathbf{X}^{T}\mathbf{X}$ は行列で、 $\lambda$ はスカラなので普通にくくることができない）

さて、式(2)では二乗和誤差に $\displaystyle\frac{\lambda}{2}\| \mathbf{w}\|^2$ を追加しましたが、正則化項は一般的に $\displaystyle\frac{\lambda}{2} \sum_{j=1}^{M}|w_j|^q$ と書け、 $q=2$ のときが式(2)に相当します。*1　そして正則化項を加えることは、係数 $w_j$ を

$\displaystyle \sum_{j=1}^{M}|w_j|^q \leq \eta \tag{4}$

の範囲に制限していることを意味します。なぜそう言えるのか？はこちらに書きました→正則化項（罰金項）の意味

さて、式(2)を使って過学習が防げるかどうか確認しました。線形回帰を最小二乗法で解くで確認した過学習の状態に対して正則化項を追加してみます。 f:id:opabinia2:20180124224859p:plain `

緑色の線が正則最小二乗法です。滑らかな近似線になって、過学習がちゃんと防げていることが確認できました！*2

今回のコードです。

import numpy as np
import matplotlib.pyplot as plt


# y = w0*x^0+....wM*x^M　を、引数xの配列数分求める
def y(w, x, M):
    X = np.empty((M + 1, x.size))
    for i in range(M + 1):
        X[i, :] = x ** i
    return np.dot(w.T, X)


# ランダムシードを固定
np.random.seed(0)
# 多項式の最大べき乗数(x^0+...+x^M)
M = 9
# 訓練データ数
N = 10

# 正則化係数λ（参考書に倣った値）
lam = np.exp(-18)

# 訓練データの列ベクトル
x = np.linspace(0, 1, N).reshape(N, 1)
# 訓練データtの列ベクトル
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)
t = t.reshape(N, 1)

# 行列Phiを作成
Phi = np.empty((N, M+1))
for i in range(M+1):
    Phi[:, i] = x.reshape(1, N) ** i

# 係数wの列ベクトルを解析的に求める
w = np.linalg.solve(np.dot(Phi.T, Phi), np.dot(Phi.T, t))

# 正則化項付きで求める
w2 = np.linalg.solve(np.dot(Phi.T, Phi) + lam*np.eye(M+1), np.dot(Phi.T, t))

# 求めた係数wを元に、新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 1, 100)
y1 = y(w, x2, M)
y2 = y(w2, x2, M)

# 結果の表示
plt.xlim(0.0, 1.0)
plt.ylim(-1.5, 1.5)
plt.scatter(x, t)
plt.plot(x2, y1.T)
plt.plot(x2, y2.T)
plt.show()

*1:ベクトル $\mathbf w$ の各要素の二乗和はノルムの二乗に等しいので。参考：ベクトルの内積と二乗和、ベクトルのノルム

*2:最初、sin(x) (-3<x<3)の曲線で試したら、正則化係数をいろいろ変えてもうまく過学習が防げませんでした。ですので手元の参考書に倣ってxの範囲を0<x<1としました。係数の値が大きくなることを防ぐ、という目的に対して、sin(x) (-3<x<3)だと過学習を起こしている時の係数もそれほど大きくなかったことが原因のように思いましたが、じゃあそういうときはどうすればいいのか？詳しい方、教えてください、、、→後日追記：λの調整でなんとかなりました。exp(-18)とは全然オーダーが違った。