Ridge回帰とLasso

Ridge回帰とLassoの簡単な方法の紹介です（メモ）．

Rのパッケージglmnetを使って簡単にできます*1．

線形回帰モデルとは
\[
y=X\beta+\varepsilon,\quad y\in\mathbf{R}^n,\quad X\in\mathbf{R}^{n\times p},\quad \beta\in\mathbf{R}^p,\quad \varepsilon\in\mathbf{R}^n,
\]
というモデルです．$y$は目的変数，$X$は説明変数（計画行列）です．誤差項$\varepsilon$には平均0の正規分布を仮定するのが基本的です．

さてデータ$(X,y)$からパラメータ$\beta$を推定することが線形回帰モデルの目標ですが，Ridge回帰は
\[
\hat{\beta}_{\mathrm{ridge}}:=\mathrm{argmin}_{\beta} ~\{(y-X\beta)^2 +\lambda |\beta|^2\},
\]
を推定量として得る方法です．ここで$|\cdot |$はユークリッドノルム，$\lambda\geq 0$は正則化項（自分で決めるもの）です．$\lambda=0$のときは最小二乗推定量となります．$\lambda>0$のとき$\beta$の大きさに罰則をつけることで過学習を防ぎます．特に説明変数の間に相関が強い場合（共線性）有効な手法です．共線性があるときに二乗ノルムの罰則を与えればよさそうという簡単な理由は，最小二乗推定量(OLS)の分散が$X^{\top}X$の行列式に比例するからです:
\[
\mathrm{E}_{\beta}|\hat{\beta}_{\mathrm{OLS}}|^2\propto \det(X^{\top}X)+|\beta|^2
\]
Hoerl and Kennard (1970)によって提案されたよく知られた手法です．

Lassoは
\[
\mathrm{argmin}_{\beta} ~\{(y-X\beta)^2 +\lambda |\beta|\},
\]
を推定量として得る方法です．Ridgeと違い$L^1$の罰則をつけることにより疎（スパース）な解を得ることができます．つまりLassoで得られた推定量$\hat{\beta}$の成分の多くが真に0をとりやすい，という性質があります．Tibshirani(1996)により提案され，非常によく研究されている手法の一つです．

どちらの手法も最小二乗法に比べて良い推定や予測ができるということがあります．しかし実際にridgeとlassoをやるには正則項$\lambda$を決める必要があります．正則化項の基本的な決め方の方法に，交差確認による決め方があります．

まず$\lambda$の候補の値（複数個，たくさん）を用意します．学習データを$n$個に分割し，そのうちの$n-1$個のパートと候補の$\lambda$の値を基に$\beta$を推定します($\hat{\beta}$)．残しておいた分割に対して，学習した$\hat{\beta}$を用いた予測値をだし，実際の値と予測値との誤差を計算します．これを$n$回繰り返し，一番誤差が小さかった$\lambda$を選ぶ方法です．手間がかかります．

実際あるデータでやってみた結果です：

f:id:mkprob:20141224000310p:plain

横軸は候補の$\lambda$の値（のlog）を表し，縦軸はその$\lambda$の下での予測誤差をプロットしたものです．
赤線は予測誤差の平均，それと上下のバウンドがついています．
縦の値が一番小さい$\lambda$を選択します．この図だと2.225(logだと0.799)となっています．

実データでの例

RのパッケージSMPracticalsのairpollution dataを使った実験をします．
このデータは$n=60$，$p=15$のデータです．変数の数を増やすために二次の項（変数$i$と変数$j$の積）を加えて$p=135$としたデータでRidgeとLassoを比較します．100回ランダムに学習データとテストデータを半分ずつに分けて予測誤差を比較しました．$\lambda$を決める際の分割数はデフォルトの10としています．

#Lasso and Ridge
set.seed(1) #random partition

#input data
airpollution <- as.matrix(read.table("air_pollution.txt"))
airp.X <- airpollution[,-16]
airp.Y <- airpollution[,16]

#high dimensional
#二次の項（X_i*X_j）の追加
airp.X.low <- airp.X
airp.X <- matrix(0,ncol=135,nrow=nrow(airp.X))
airp.X[,1:15] <- airp.X.low
airp.X[,16:30] <- airp.X.low^2
for(i in 1:105){
	airp.X[,30+i] <- apply(airp.X[,combn(15,2)[,i]],1,prod)
}

#grid: 正則化の探す範囲（グリッドサーチ）
grid <- 10^seq(8,-3,length=1000)
#MSE
num.crossval <- 100
mse.mat <- matrix(0,nrow=num.crossval,ncol=2)

for(i in 1:num.crossval){
	train <- sample(nrow(airp.X),30) #ランダムに学習データ（30個）とテストデータ（30個）に分割
	test <- -train

	#ridge
	cv.out.ridge <- cv.glmnet(airp.X[train,],airp.Y[train],alpha=0,lambda=grid) #交差確認
	bestlam.ridge <- cv.out.ridge$lambda.min #最適な正則化項の値
	ridge.mod <- glmnet(airp.X[train,],airp.Y[train],alpha=0,lambda=bestlam.ridge)
	ridge.pred <- predict(ridge.mod,s=bestlam.ridge,newx=airp.X[test,]) #予測
	mse.mat[i,1] <- sqrt(mean((ridge.pred-airp.Y[test])^2))

	#lasso
	cv.out.lasso <- cv.glmnet(airp.X[train,],airp.Y[train],alpha=1,lambda=grid) #交差確認
	bestlam.lasso <- cv.out.lasso$lambda.min #最適な正則化項の値
	lasso.mod <- glmnet(airp.X[train,],airp.Y[train],alpha=1,lambda=bestlam.lasso)
	lasso.pred <- predict(lasso.mod,s=bestlam.lasso,newx=airp.X[test,]) #予測
	mse.mat[i,2] <- sqrt(mean((lasso.pred-airp.Y[test])^2))
}

#plot
boxplot(mse.mat,names=c("Ridge","Lasso"))