Ridge回帰の使い所を考える

はじめに

本記事では回帰係数の推定方法としてのRidge回帰（L2正則化）の使い所について考えたいと思います。
結論から言うと、サンプルサイズが不十分な状況下でRidge回帰を用いてより真の値に近い回帰係数を得る確率を高めるという使い方ができそうです。

OLS推定量とRidge推定量

簡単のために、以下の線形回帰モデルを考えます。

\begin{align}
Y = X\beta + \epsilon
\end{align}

ここで、 $X$ は $N \times K$ の行列で、 $\beta$ は $K\times 1$ の列ベクトルとします。 $N$ はサンプルサイズ、 $K$ は特徴量の数になります。
このとき、OLS推定量は二乗誤差を最小にするような $b$ で与えられます。
\begin{align}
\widehat{\beta}_{\mathrm{OLS}} &= \arg\min_b \left(Y - X b\right)^\top\left(Y - X b\right) \\
&= \left(X^\top X \right)^{-1}X^\top Y
\end{align}

標準的な仮定のもとで、OLS推定量は推定量 $\widehat{\beta}_{\mathrm{OLS}}$ は真の値 $\beta$ の一致推定量となることが知られています。*1

これに対して、 $b$ のとれる値に制限をかけることで過学習を抑えようとするのがRidge回帰になります。
具体的にはOLSの目的関数に正則化項 $\lambda||b||^2$ を付け加えたものを最小化します。 $\lambda$ は正則化の強さを表すハイパーパラメータになります。

\begin{align}
\widehat{\beta}_{\mathrm{Ridge}} &= \arg\min_b \left(Y - X b\right)^\top\left(Y - X b\right) + \lambda||b||^2 \\
&= \left(X^\top X + \lambda I_K \right)^{-1}X^\top Y
\end{align}

ここで $I_K$ は $K\times K$ の単位行列です。
$\widehat{\beta}_{\mathrm{OLS}}$ と $\widehat{\beta}_{\mathrm{Ridge}}$ を比較すれば明らかですが、両者の違いは $\lambda I_K$ のみです。
なので、 $\lambda = 0$ のときOLS推定量とRidge推定量は完全に一致します。反対に $\lambda$ を大きくするとRidge推定量は0方向に引っ張られます。
つまり、Ridge回帰はパラメータが極端な値を取らないよう安定させる（過学習を抑える）というリターンを得るために、（多少の）バイアスがかかるという犠牲を払うモデルであると言えます。
これは機械学習の文脈でバイアスとバリアンスのトレードオフとしてよく知られています。

シミュレーション

それでは実際にOLSとRidge回帰の挙動の違いについてシミュレーションで確認してみましょう。
まずは必要なパッケージを読み込みます。*2

import pandas as pd
import numpy as np

from sklearn.linear_model import LinearRegression, Ridge,  RidgeCV
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_spd_matrix #分散共分散行列を作る

from plotnine import * # 可視化用。

ちょっと極端な方が結果がクリアに見えると思うので、サンプルサイズは100しかないんだけど80個の特徴量について回帰係数を推定したい、というシチュエーションを想定することにします。
下記のコードではランダムにデータを生成してOLSとRidge回帰で回帰係数を推定する、ということを繰り返しています。

np.random.seed(42)
I = 1000
N = 100
K = 80

beta = np.random.uniform(-1, 1, K)

betas = np.zeros((I, 2, K))
for i in range(I):
    X = np.random.multivariate_normal(np.zeros(K), make_spd_matrix(K), N)
    X = StandardScaler().fit_transform(X)
    U = np.random.normal(0, 1, N)
    Y = X @ beta + U
    
    if i == 0: #一回目のループでRidgeのハイパーパラメータを決める
        ridge_alpha = RidgeCV(alphas=np.random.uniform(0.1, 20, 100), cv=5).fit(X, Y).alpha_
        print(f'ridge: {ridge_alpha}')
    
    beta_ols = LinearRegression().fit(X, Y).coef_
    
    beta_ridge = Ridge(alpha=ridge_alpha).fit(X, Y).coef_
    
    betas[i, :, :] = np.concatenate([beta_ols[:, np.newaxis], beta_ridge[:, np.newaxis]], axis=1).T
    
    if i % 200 == 0:
        print(i)

# 可視化のためにシミュレーションの結果をデータフレームとして持ち直す
df = None
for i in range(I):
    tmp = (pd.DataFrame(data = betas[i, :, :].T, columns=['OLS', 'Ridge'])
           .assign(Truth = beta, k = lambda df: df.index, i = i))
    
    if df is None:
        df = tmp
    else:
        df = df.append(tmp)
        
df = df.melt(id_vars=['i', 'k', 'Truth'], value_vars=['OLS', 'Ridge'], value_name='coefficient', var_name='method')

シミュレーションの結果を可視化してみましょう。

g = (ggplot(aes('factor(k)', 'coefficient', fill='method', color='method'), df[df.k.between(1, 10)]) +
     geom_hline(yintercept=0, size=1, color='gray') +
     geom_violin(alpha=0.5) +
     geom_point(aes('factor(k)', 'Truth'), size=4, color='black', fill='black') +
     scale_fill_brewer(name='Method', type='qual',palette='Dark2') +
     scale_color_brewer(name = 'Method', type='qual',palette='Dark2') + 
     lims(y=(-3, 3)) + 
     labs(x='beta_hat(k)', y='Coefficient') + 
     theme_bw() +
     theme(figure_size=(15, 7)))

print(g)

f:id:dropout009:20190102204239p:plain

図は初めの10個の回帰係数を取り出して、シミュレーション1000回分の回帰係数の分布を確認しています。
縦軸が回帰係数の値、横軸が回帰係数の番号です。
真ん中の黒い点が真の回帰係数の値で、緑がOLSで求めた回帰係数の分布、オレンジがRidge回帰での分布になります。

図からは各手法で求めた回帰係数について、以下の性質が見て取れます。

OLSの場合はバイアスがなく真の値を中心に分布しているが、その分散は大きく、見当はずれの値を取ることもある
Ridge回帰の場合はゼロ方向にバイアスがかかっている。その一方で回帰係数は非常に安定しており、真の値に近い部分に密度高く分布しており、極端な値を取ることもない。

Ridge回帰の使いどころ

数式とシミュレーションの結果から、サンプルサイズが不十分な状況下でRidge回帰を用いてより真の値に近い回帰係数を得る確率を高めるという使い方が考えられそうです。*3
特に、特徴量のインパクトを見える際により保守的な推定をしたい場合には有効と言えそうです。Ridge回帰は極端な回帰係数をとらず、またRidge回帰のバイアスは0方向にかかるので、Ridge回帰ではより控え目なインパクトが推定されると言えるかと思います。

*1:https://www.ssc.wisc.edu/~bhansen/econometrics/

*2:ちなみにplotnineはRの代表的な可視化パッケージであるggplot2と同じ記法が使えるpython用のパッケージです。若干制限もありますが、ggplot2に慣れたRユーザーにとっては使いやすいかと思います。

*3:サンプルサイズが十分な状況ではOLSで一致推定を行うのが妥当かと思います。