過学習とL2正則化

## 過学習

過学習とは、モデルが学習データにのみフィットしており、未知のデータの予測精度が低くなること。一般にモデルを複雑にするほど過学習は起きやすい。
未知のデータを予測する能力を「汎化能力」と言い、未知のデータを予測した際の誤りを「汎化誤差」と言う。

### １．１．未知のデータを作る。
まず、未知のデータを作って過学習を実際に見てみる。

```py
import numpy as np
import math
import matplotlib.pyplot as plt

np.random.seed(seed=32)
# create samples
sample_size = 50
test_size = 20
noise_size = 0.2

dataX = np.linspace(0.0, 5.0, num=1000).reshape((1000, 1))
# create train data
x = np.random.permutation(dataX)[:sample_size]
noise = noise_size * np.random.randn(sample_size, 1)
func = np.vectorize(math.sin)
y = func(x * 1.6) + noise
plt.plot(x, y, 'o', label='train data')

# create test data
testx = np.random.permutation(dataX)[test_size * -1:]
noise = noise_size * np.random.randn(test_size, 1)
testy = func(testx * 1.4) + noise
plt.plot(testx, testy, 'ro', label='test data')
```

全データの先頭50件を訓練データとし、末尾20件をテストデータとしている。
*test_yを作成するとき、正弦関数の係数をわざと小さめに調整している。
（実際は$y$が異なる分布となるので同一集団ではなくなるが、ここでは気にせず先に進める。）


これに前回の学習モデルを適用した結果と一緒に図示すると、
<img width="50%" src="https://qiita-image-store.s3.amazonaws.com/0/204732/52bd639e-b14f-b517-ff9c-0e29b7998410.png">

汎化能力に関しては、図からどっちがどうか分かりにくい。

### １．２．汎化誤差を計算する。
```py
dim = [3, 9]
for d in dim:
    # make predictions
    poly = PolynomialFeatures(degree=d)
    X = poly.fit_transform(x)
    testX = poly.fit_transform(testx)
    XTX = np.dot(X.T, X)
    XTX_inv = np.linalg.inv(XTX)
    a = np.dot(XTX_inv, np.dot(X.T, y))
    print "bias", np.sum((y - np.dot(X, a)) ** 2) / sample_size
    print "generalization", np.sum((testy - np.dot(testX, a)) ** 2) / test_size
    Xt = poly.fit_transform(dataX)
    yt = np.dot(Xt, a)
    plt.plot(dataX, yt, label='%d polynomial' % d)
    >> bias 0.163705879935
    >> generalization 0.28673629963
    >> bias 0.0444546042927
    >> generalization 0.329093164686
```
バイアスは９次元の方が小さいが、汎化誤差は３次元の方が小さくなっている。

## L2（リッジ）正則化

次に過学習を抑える手法である正則化を取り入れる。
損失関数にパラメータの２乗ノルムを加える。

```math
 E(a)= \parallel y-Xa \parallel ^{2}+\lambda\parallel a \parallel ^{2}
```
正則化項により、aの値が大きくなるのを防ぐ。
（大きくなると、損失関数が大きくなるため）
＊$\lambda$は正則化パラメータ（大きくするほど正則化の度合いも高くなる。）

損失関数を最小化する$a$を求めるため、$E$を$a$で微分すると、

```math
a=(X^TX+\lambda I)^{-1}X^Ty
```
となる。

### ２．１．正則化の効果を確認する。

正則化パラメータに「**10, 1, 0, 0.1**」を用意して効果を確認する。（0は正則化なし）
<img width="50%" src="https://qiita-image-store.s3.amazonaws.com/0/204732/dffbaecb-654d-89c7-a5c6-dbf66edd3e6a.png">


図の左側を見ると違いが分かりやすい。
正則化パラメータを**10**にした場合が、最もモデルの線形が滑らかになっている。


```py
d = 9
L2params = [10, 1, 0, 0.1]
for p in L2params:
    # make predictions
    poly = PolynomialFeatures(degree=d)
    X = poly.fit_transform(x)
    testX = poly.fit_transform(testx)
    XTX = np.dot(X.T, X)
    L2im = np.identity(d + 1)*p
    L2im[0,0] = 0
    XTX = XTX + L2im
    XTX_inv = np.linalg.inv(XTX)
    a = np.dot(XTX_inv, np.dot(X.T, y))
    print ">> lambda %.1f :" % p, np.sum(a**2)
    print ">> bias :", np.sum((y - np.dot(X, a)) ** 2) / sample_size
    print ">> generalization :", np.sum((testy - np.dot(testX, a)) ** 2) / test_size
    Xt = poly.fit_transform(dataX)
    yt = np.dot(Xt, a)
    plt.plot(dataX, yt, label='regularization rate %.1f' % p)

>> lambda 10.0 : 0.433173909089
>> bias : 0.0813137786099
>> generalization : 0.322284436323
>> lambda 1.0 : 0.588098400133
>> bias : 0.0599131733257
>> generalization : 0.308749806939
>> lambda 0.0 : 290.780236508
>> bias : 0.0444545966264
>> generalization : 0.32903307304
>> lambda 0.1 : 2.07354251377
>> bias : 0.050158588753
>> generalization : 0.310930705156

```
また、正則化なしの場合がパラメータのL2ノルムが最も大きい(290)ことが分かる。
但し、今回は汎化性能は僅かな改善しか見られなかった。。

### ２．２．正則化の効果を考える。
正則化により、パラメータがより滑らかになることが分かった。
その理由についてもう少し考えてみる。

こんなデータセットでの回帰を考える。

| No |  説明変数1 | 説明変数2 | 目的変数 |
|:-:|:-:|:-:|:-:|
| 1  | 3  | 5  | 5  |
| 2  | 4  | 2  | 3  |
| 3  | 3  | 7  | 8  |
| 4  | 4  | 8  | 9 |
これを行列で表すと、

```math
X=
\begin{pmatrix}
1 & 3 & 5 \\
1 & 4 & 2 \\
1 & 3 & 7 \\
1 & 4 & 8 
\end{pmatrix}
\;y=
\begin{pmatrix}
5 \\
3 \\
8 \\
9  
\end{pmatrix}
```
となる。
$(X^TX)^{-1}X^Ty$を一つずつ確認する。
まず、$X^TX$であるが、これは計算すると、

```math
X^TX=
\begin{pmatrix}
4 & 14 & 22 \\
14 & 50 & 76 \\
22 & 76 & 142 
\end{pmatrix}
```
となる。
逆行列$(X^TX)^{-1}$は、

```math
(X^TX)^{-1}=
\begin{pmatrix}
16.55 & -3.95 & -0.45 \\
-3.95 & 1.05 & 0.05 \\
-0.45 & 0.05 & 0.05 
\end{pmatrix}
```
となる。
$a$を求める際、$(X^TX)^{-1}$の各行が$a$の各変数に影響する。（１行目は$a$の１つ目の変数に影響する）
$X^Ty$は正則化による影響を受けないので計算は割愛する。

ここで、正則化項（$\lambda=1$）を導入して計算にどう影響しているかを確認する。

```math
\lambda I=
\begin{pmatrix}
0 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 
\end{pmatrix}
```
（１行目は説明変数に関係ないパラメータなので含めない）

```math
X^TX+\lambda I=
\begin{pmatrix}
4 & 14 & 22 \\
14 & 51 & 76 \\
22 & 76 & 143 
\end{pmatrix}
```
これだけ見ると、正則化の効果はほとんど無いように感じる。
しかし、$(X^TX+\lambda I)^{-1}$の結果は、

```math
(X^TX+\lambda I)^{-1}=
\begin{pmatrix}
8.819 & -1.918 & -0.337 \\
-1.918 & 0.511 & 0.023 \\
-0.337 & 0.023 & 0.988 
\end{pmatrix}
```
と変化しており、全体的に値が小さくなっていることが分かる。
これが正則化による効果であり、$a$のL2ノルムが小さくなる理由である。

### ２．３．scikit-learnで確認する。
最後に「linear_model」ライブラリとの一致を確認する。
scikit-learnの場合、

```py
from sklearn import linear_model
crf = linear_model.Ridge(alpha=p)
crf.fit(X, y)
```
でパラメータが算出可能。

```py
for p in L2params:
    # make predictions
    poly = PolynomialFeatures(degree=d)
    X = poly.fit_transform(x)
    testX = poly.fit_transform(testx)
    XTX = np.dot(X.T, X)
    L2im = np.identity(d + 1) * p
    L2im[0, 0] = 0
    XTX = XTX + L2im
    XTX_inv = np.linalg.inv(XTX)
    a = np.dot(XTX_inv, np.dot(X.T, y))
    print "compute by myself ------"
    print("intercept : ", a[0])
    print("coef : ", a[1:].T)
    crf = linear_model.Ridge(alpha=p)
    crf.fit(X, y)
    print "compute by sklearn ------"
    print("intercept : ", crf.intercept_)
    print("coef : ", crf.coef_)
plt.legend()
compute by myself ------
('intercept : ', array([ 0.641]))
('coef : ', array([[ 0.101,  0.063, -0.007, -0.064, -0.051,  0.033, -0.002, -0.001,  0.   ]]))
compute by sklearn ------
('intercept : ', array([ 0.641]))
('coef : ', array([[ 0.   ,  0.101,  0.063, -0.007, -0.064, -0.051,  0.033, -0.002,
        -0.001,  0.   ]]))
compute by myself ------
('intercept : ', array([ 0.394]))
('coef : ', array([[ 0.503,  0.27 , -0.057, -0.24 , -0.11 ,  0.172, -0.06 ,  0.009, -0.   ]]))
compute by sklearn ------
('intercept : ', array([ 0.394]))
('coef : ', array([[ 0.   ,  0.503,  0.27 , -0.057, -0.24 , -0.11 ,  0.172, -0.06 ,
         0.009, -0.   ]]))
compute by myself ------
('intercept : ', array([ 0.219]))
('coef : ', array([[ -1.427,   9.773, -12.472,   6.066,  -0.835,  -0.344,   0.159,
         -0.024,   0.001]]))
compute by sklearn ------
('intercept : ', array([ 0.219]))
('coef : ', array([[  0.   ,  -1.427,   9.772, -12.472,   6.065,  -0.834,  -0.344,
          0.159,  -0.024,   0.001]]))
compute by myself ------
('intercept : ', array([ 0.099]))
('coef : ', array([[ 1.228,  0.316, -0.49 , -0.382,  0.262, -0.014, -0.018,  0.004, -0.   ]]))
compute by sklearn ------
('intercept : ', array([ 0.099]))
('coef : ', array([[ 0.   ,  1.228,  0.316, -0.49 , -0.382,  0.262, -0.014, -0.018,
         0.004, -0.   ]]))

```

scikit-learnは、coefの最初に0が入っているが、多分$X$の１列目に１がある影響だと思う。
それ以外は一致しているのでOKとする。