六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

RでL1 / L2正則化を実践する

L1 / L2正則化と言えば機械学習まわりでは常識で、どんな本を見てもその数式による表現ぐらいは必ず載ってる*1わけですが、そう言えばあまり実務では真面目にL1 / L2正則化入れてないなと思ったのと、Rでやるなら普通どうするんだろう?と思ったので*2、もはや周回遅れみたいな感じではありますが備忘録的に実践してみようかと思います。

L1 / L2正則化って何だっけ


ということで復習(自分の記憶が合っているかどうかの確認)。。。PRMLにも載ってる有名な図がありますが、あれの説明が直感的には最も分かりやすいかと思います。これは重み付けベクトルが2次元の場合、つまりw_1w_2を求めるという問題を想定した図です。もうちょっと言えば2次元データに対する分類or回帰問題ということで。


f:id:TJO:20150228165325p:plain


基本的には分類器にせよ回帰モデルにせよ、学習データに対する誤差*3E_D(w)を定式化して、これを最小化するように(このケースでは2つの)重み付けを最適化していくという問題を解いています。なのですが、闇雲に誤差ばかり最小化すると多くの機械学習の本にも書いてあるように過学習(overfitting)に陥る恐れがあります。


そこでこれにペナルティ項E(w)を加えて、E_D(w) + \lambda E(w)を最小化するように最適化問題を解く問題に置き換えることで、「ほどよく」過学習を避けつつ穏当なモデルに落ち付かせようというのが、いわゆる「正則化」(regularization / penalized regression)の考え方です。


上の図はその様子を模式化したもので、要は w_1, w_2からなる重み付けベクトルの空間上に学習データに対する誤差項の広がりを配置しこれを動かすことで誤差を最小化させるわけですが、そこにペナルティ項を付けておくことで w_1, w_2の定まり方が変わるというお話です。左側のL1ノルム正則化であれば尖った形の領域に接するようにしなければならないため自ずとw_2の軸上に最適点が定まってw_1 = 0となり、右側のL2ノルム正則化であれば円状の領域に接するようにしなければならないためw_1, w_2ともにnon-zeroな値が入る、という感じになるわけです。


なお『イラストで学ぶ 機械学習 最小二乗法による識別モデル学習を中心に (KS情報科学専門書)』(通称杉山本)pp.31-35及びpp.41-47にある通り、L1正則化項は||\mathbf{\theta}||_1 \leq R, ||\mathbf{\theta}||_1 = \displaystyle \sum^b_{j=1} | \theta_j|のように重み付けベクトル\mathbf{\theta}の要素の絶対値の総和がR以下という領域、すなわちカクカクした四角形のような領域で表されます。またL2正則化項は||\mathbf{\theta}||^2 \leq Rのように重み付けベクトル\mathbf{\theta}の要素の二乗和がR以下という領域、すなわち円形の領域で表されます。


で、L1正則化の場合はこの図の例のように「w_2だけに値が入りw_1はゼロになる」というような、いわゆるスパース(疎)な解を得やすいことが分かります。そこでL1正則化「不要なパラメータを削りたい」(次元・特徴量削減)という時によく使われるんですね。一方でL2正則化の場合は過学習を抑えて汎化された滑らかなモデルを得やすいことが知られています。なおL1正則化回帰はLasso回帰、L2正則化回帰はRidge回帰とも呼ばれ、教科書によっては主にこちらの名前で書かれているものもあります*4


ちなみにL1 / L2正則化項は数学的には単なる線形和で表せるので、普通にL1を0.6, L2を0.4というように割合を定めてミックスさせることもできます。これをElastic net正則化と呼びます。


(※間違ってたらいつも通り炎上ラーニング大歓迎ですのでビシバシ突っ込んで下さい)


RでL1 / L2正則化回帰するなら{glmnet}で


ということでRにはこの辺全部実装してあるナイスなパッケージ{glmnet}があるので、これを使って色々試してみようかと思います。なお{glmnet}にはご丁寧にwebページ版のvignetteが用意されているので、これを見ながらやればOKかと。


Glmnet Vignette - Trevor Hastie and Junyang Qian, Stanford, June 26, 2014


データセットぐらいは違うものを使いたい気もするので、ここでは正規線形モデル向けにUCI Machine Learning Repositoryから"Energy Efficiency"データセットを、ロジスティック回帰向けに同じくテニス四大大会スコア&パフォーマンスデータセットを使うことにします。


ということで以下のGitHubリポジトリからそれぞれ "energey_efficency.csv" と "men.txt" を落としてきていただきたく。それぞれd1,d2という名前のデータフレームにして読み込んでおきましょう。



まずは{glmnet}を入れます。インストールしていない人はインストールもしておきましょう。

> install.packages("glmnet")
> library(glmnet)


これで準備が整いました。では実際のデータに適用してみましょう。


Energy Efficiencyデータセットの場合


まず"Energy Efficiency"データセットd1について。これは様々な住宅における暖房・冷房効率の良し悪しを、個々の住宅の「寸法」「面積」などによってモデリングしようというものです(9列目が暖房効率、10列目が冷房効率なので分けてモデリングする必要がある)。


ということで、特に何も考えず以下のようにやってみます。最初にlm関数で普通の線形回帰をやってから、その次にglmnet関数でL1正則化Lasso回帰をかけます。

> d1.lm1<-lm(heating_load~.,d1[,-10])
> summary(d1.lm1)

Call:
lm(formula = heating_load ~ ., data = d1[, -10])

Residuals:
    Min      1Q  Median      3Q     Max 
-9.8965 -1.3196 -0.0252  1.3532  7.7052 

Coefficients: (1 not defined because of singularities)
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                84.014521  19.033607   4.414 1.16e-05 ***
relative_compactness      -64.773991  10.289445  -6.295 5.19e-10 ***
surface_area               -0.087290   0.017075  -5.112 4.04e-07 ***
wall_area                   0.060813   0.006648   9.148  < 2e-16 ***
roof_area                         NA         NA      NA       NA    
overall_height              4.169939   0.337990  12.337  < 2e-16 ***
orientation                -0.023328   0.094705  -0.246  0.80550    
glazing_area               19.932680   0.813986  24.488  < 2e-16 ***
glazing_area_distribution   0.203772   0.069918   2.914  0.00367 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.934 on 760 degrees of freedom
Multiple R-squared:  0.9162,	Adjusted R-squared:  0.9154 
F-statistic:  1187 on 7 and 760 DF,  p-value: < 2.2e-16

> d1.glmnet1<-glmnet(as.matrix(d1[,-c(9,10)]),as.matrix(d1[,9]),family="gaussian",alpha=1)
> plot(d1.glmnet1)
> coef(d1.glmnet1)
9 x 73 sparse Matrix of class "dgCMatrix"
   [[ suppressing 73 column names ‘s0’, ‘s1’, ‘s2’ ... ]]
                                                                           
(Intercept)               22.3072 19.9169408 17.7390251 15.754589 13.946446
relative_compactness       .       .          .          .         .       
surface_area               .       .          .          .         .       
wall_area                  .       .          .          .         .       
roof_area                  .       .          .          .         .       
overall_height             .       0.4552876  0.8701286  1.248116  1.592525
orientation                .       .          .          .         .       
glazing_area               .       .          .          .         .       
glazing_area_distribution  .       .          .          .         .       
                                                                                 
(Intercept)               12.298932 10.797780 9.429985 8.183702 7.048135 6.013448
relative_compactness       .         .        .        .        .        .       
surface_area               .         .        .        .        .        .       
wall_area                  .         .        .        .        .        .       
# 以下略

f:id:TJO:20150228174004p:plain


glmnet関数は説明変数と目的変数を分けてそれぞれmatrix型で与えなければならない点に要注意。family引数に"gaussian"を指定すると線形回帰してくれます。alphaは1ならL1正則化(Lasso)、0ならL2正則化(Ridge)、その間の値ならelastic netになります。


で、結果を見ると。。。plotして見られるのがL1ノルムの値に依存した個々の偏回帰係数の変化の様子で、coefで見られるのがL1ノルムの値ごとに並べられた偏回帰係数のリストです。でもこれだと色々あり過ぎてよく分からないですよね。


基本的に{glmnet}で扱うL1 / L2正則化回帰は、L1 / L2ノルムの値に依存して偏回帰係数の推定結果が変わります。なのでパラメータチューニングが必要なのですが、{glmnet}にはcross validationを行ってチューニングをサポートしてくれる関数があるので、やってみようと思います。

> d1.cv.glmnet1<-cv.glmnet(as.matrix(d1[,-c(9,10)]),as.matrix(d1[,9]),family="gaussian",alpha=1)
# heating_loadのLasso回帰
> d1.cv.glmnet2<-cv.glmnet(as.matrix(d1[,-c(9,10)]),as.matrix(d1[,10]),family="gaussian",alpha=1)
# cooling_loadのLasso回帰
> plot(d1.cv.glmnet1)
> plot(d1.cv.glmnet2)
> coef(d1.cv.glmnet1,s=d1.cv.glmnet1$lambda.min)
9 x 1 sparse Matrix of class "dgCMatrix"
                                     1
(Intercept)                44.22794671
relative_compactness      -42.97744582
surface_area                .         
wall_area                   .         
roof_area                  -0.10509037
overall_height              4.63150008
orientation                -0.01343843
glazing_area               19.86421000
glazing_area_distribution   0.19789049
> coef(d1.cv.glmnet2,s=d1.cv.glmnet2$lambda.min)
9 x 1 sparse Matrix of class "dgCMatrix"
                                     1
(Intercept)                82.60593418
relative_compactness      -62.81477092
surface_area               -0.03388986
wall_area                   .         
roof_area                  -0.08296767
overall_height              4.46183969
orientation                 0.11939088
glazing_area               14.70239391
glazing_area_distribution   0.03943679

f:id:TJO:20150228175130p:plain

f:id:TJO:20150228175138p:plain


結果として、surface_areaやwall_areaが不要という結論に達することが分かりました。こんな感じでL1 / L2正則化回帰を{glmnet}では実践することができます。ただしelastic net時のLasso : Ridgeの割合を決めるためのチューニング手法は実装されていないようなので、そこは色々考える必要がありそうです。。。


テニス四大大会データセットの場合


これは以前の記事を参考にしていただければと。ロジスティック回帰であっても以下のようにfamily引数を"binomial"にすればやれます。ついでに"women.txt"をデータフレームdwとして読み込んでおいて、以下のようにテスト分類もやってみましょう。とりあえずここでは以前の記事同様に次元削減を念頭に置いているため、L1正則化に限ることとします。

> d2<-d2[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]
> d2.cv.glmnet<-cv.glmnet(as.matrix(d2[,-1]),as.matrix(d2[,1]),family="binomial",alpha=1)
> plot(d2.cv.glmnet)
> coef(d2.cv.glmnet,s=d2.cv.glmnet$lambda.min)
25 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  5.535541e-02
FSP.1        5.252102e-02
FSW.1        1.287456e-01
SSP.1       -3.918233e-05
SSW.1        1.703759e-01
ACE.1        .           
DBF.1       -9.703481e-02
WNR.1        4.140972e-02
UFE.1       -8.009685e-03
BPC.1        3.722087e-01
BPW.1        2.230594e-01
NPA.1        .           
NPW.1        .           
FSP.2       -3.679382e-02
FSW.2       -1.733349e-01
SSP.2        .           
SSW.2       -1.513321e-01
ACE.2        1.601675e-02
DBF.2        5.009638e-02
WNR.2       -2.295361e-02
UFE.2       -2.670256e-03
BPC.2       -3.579647e-01
BPW.2       -2.164438e-01
NPA.2        .           
NPW.2        1.463902e-02

> dw<-dw[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]
> table(dw$Result,round(predict(d2.cv.glmnet,as.matrix(dw[,-1]),s=d2.cv.glmnet$lambda.min,type='response'),0))
   
      0   1
  0 216  11
  1  19 206

> sum(diag(table(dw$Result,round(predict(d2.cv.glmnet,as.matrix(dw[,-1]),s=d2.cv.glmnet$lambda.min,type='response'),0))))/nrow(dw)
[1] 0.9336283
# お、正答率が前回の記事超えてる!

f:id:TJO:20150228180721p:plain


CVプロットが思ったよりはっきりとした極小値を見せているので、これはL1正則化の効果があるだろうな~と思ったらズバリ出ました。しかもテストデータとして女子データセットの分類をかけてみたら、ステップワイズ法で説明変数を絞り込んだ時よりも正答率が高いじゃないですか! これは面白いですね~。


加えて興味深いのが、stepAIC関数でステップワイズ法を用いて絞り込んだ説明変数のリストと、cv.glmnet関数でL1正則化によって絞り込んだ説明変数のリストが、想像以上に食い違っている点。この辺は何かしらの理論的背景があるんだろうなぁとも思うんですが、ここでは深入りは避けておきますよということで。


Rで正則化項を含む分類器を使いたい場合


回帰は{glmnet}で片付いたんですが、分類器の場合は{LiblineaR}で線形SVMを使うことができます。が、線形SVM使うくらいならロジスティック回帰でいいんじゃ?という気が。。。


一応、そのまんまLiblineaR関数のtype引数を以下のように設定することでL1 / L2正則化を導入できます。

type


LiblineaR can produce 10 types of (generalized) linear models, by combining several types of loss functions and regularization schemes. The regularization can be L1 or L2, and the losses can be the regular L2-loss for SVM (hinge loss), L1-loss for SVM, or the logistic loss for logistic regression. The default value for type is 0. See details below. Valid options are:


for multi-class classification

0 – L2-regularized logistic regression (primal)

1 – L2-regularized L2-loss support vector classification (dual)

2 – L2-regularized L2-loss support vector classification (primal)

3 – L2-regularized L1-loss support vector classification (dual)

4 – support vector classification by Crammer and Singer

5 – L1-regularized L2-loss support vector classification

6 – L1-regularized logistic regression

7 – L2-regularized logistic regression (dual)


for regression

11 – L2-regularized L2-loss support vector regression (primal)

12 – L2-regularized L2-loss support vector regression (dual)

13 – L2-regularized L1-loss support vector regression (dual)


ということで、type=3 or 5を試してみることにします。ここではロジスティック回帰の当てはまりが良かったテニス四大大会データセットの方にだけ当てはめてみましょう。

> install.packages("LiblineaR")
> library(LiblineaR)
> d2.svm.l2l1<-LiblineaR(as.matrix(d2[,-1]),as.vector(d2[,1]),type=3)
> d2.svm.l1l2<-LiblineaR(as.matrix(d2[,-1]),as.vector(d2[,1]),type=5)
> table(dw$Result,predict(d2.svm.l2l1,newx=as.matrix(dw[,-1]),decisionValues=T)$predictions)
   
      0   1
  0 195  32
  1  10 215

> sum(diag(table(dw$Result,predict(d2.svm.l2l1,newx=as.matrix(dw[,-1]),decisionValues=T)$predictions)))/nrow(dw)
[1] 0.9070796
> table(dw$Result,predict(d2.svm.l1l2,newx=as.matrix(dw[,-1]),decisionValues=T)$predictions)
   
      0   1
  0 213  14
  1  16 209

> sum(diag(table(dw$Result,predict(d2.svm.l1l2,newx=as.matrix(dw[,-1]),decisionValues=T)$predictions)))/nrow(dw)
[1] 0.9336283
# おお、正答率がLassoロジスティック回帰と同じになった!


ということで、確かにSVMでも*5L1正則化*6をかけることでパフォーマンスが向上することが分かりました。


おまけ


とは言え、この結果を見て「どうもこのテニス四大大会データセットは線形分離可能っぽいし、ひょっとして普通の線形SVMの方がガウシアンカーネルSVMよりもパフォーマンスが良かったりするんじゃね?」という気がしたので、試してみました。

> d2$Result<-as.factor(d2$Result)
> d2.tune<-tune.svm(Result~.,data=d2,kernel='linear')
> d2.tune$best.model

Call:
best.svm(x = Result ~ ., data = d2, kernel = "linear")


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  1 
      gamma:  0.04166667 

Number of Support Vectors:  97

> d2.lsvm<-svm(Result~.,d2,kernel='linear',cost=1,gamma=d2.tune$best.model$gamma)
> table(dw$Result,predict(d2.lsvm,newdata=dw[,-1]))
   
      0   1
  0 214  13
  1  16 209

> sum(diag(table(dw$Result,predict(d2.lsvm,newdata=dw[,-1]))))/nrow(dw)
[1] 0.9358407
# ベストパフォーマンスだ。。。


うわあああああ、やっぱり線形SVMというのが重要だったんじゃんかorzorzorz 初めから線形分類器を用意すべきだったのかもしれません。。。しかもサポートベクターの数も少なくてスパース性という点でも遥かに優れてるし。。。


ということでL1 / L2正則化もいいけど、ロジスティック回帰で高正答率が得られるような線形分離可能っぽいデータセットに対しては、ランダムフォレストはむしろ不適かも&SVM使うならカーネルの選択から考えた方がいいかもしれませんよ、というレッスンでしたとさ。


参考資料


PRMLとかだと読みづらいので、適当にweb上の資料を漁りました。以下そのリストを。

あとは杉山本ですね。

個人的にはこの辺の概念的な話の説明はPRMLやはじパタよりもこの本の方が分かりやすいと思うことが多いので*7、「イラストあまり関係ないよね」と苦笑しつつも時々有難く参照してます。

*1:PRMLもそうだけど大抵は概念図入りでLassoとRidgeの違いみたいなのが出てる

*2:h2o.deeplearningでは使ってますが

*3:経験損失

*4:サポートベクターマシン入門』がそうですね

*5:線形カーネル限定ながら

*6:ただし損失関数がL2ノルム

*7:他にも古典的な機械学習のテキストでは抜けていて現代では大事な話の多くが載っているのも特徴かと