渋谷駅前で働くデータサイエンティストのブログ

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

メモリに乗り切らない大容量データ相手にL1正則化回帰を回せる{biglasso}パッケージを試してみた

今回は小ネタ。ボサーっとCRAN Task View: Machine Learningを眺めていたらこんなものを見つけました。

これ、カラクリは簡単で単にR側にはデータのポインタしか渡さず、データの実体はストレージからオンライン(ストリーミング)で読み込むようにしているというだけなんですね。


基本的に線形回帰(gaussian)とロジスティック回帰(binomial)しかまだサポートしていませんが、後者に関しては以前論文メモ記事にしたSloresも実装されている模様です。

ということで、サクッと試してみます。


{glmnet}と比較してみる


そもそもVignetteにも同じL1正則化回帰を回せる{glmnet}との比較が出ているんですが、違うデータセットで僕もやってみました。使ったのはUCI ML repositoryのOnline News Popularityです。

> d <- read.csv('OnlineNewsPopularity/OnlineNewsPopularity.csv')
> d <- d[,-c(1,2)]
> idx <- sample(nrow(d),5000,replace=F)
> d_train <- d[-idx,]
> d_test <- d[idx,]
> X <- as.matrix(d_train[,-59])
> y <- d_train[,59]
> library(glmnet)
> t <- proc.time()
> glmfit <- cv.glmnet(X, y, family='gaussian', alpha=1)
> proc.time() - t
   ユーザ   システム       経過  
     1.821      0.402      2.233 
> library(biglasso)
> X.bm <- as.big.matrix(X)
> t <- proc.time()
> cvfit <- cv.biglasso(X.bm, y, family='gaussian', seed=71, ncores=3)
# 起動準備中&各種アラートが出る
> proc.time() - t
   ユーザ   システム       経過  
     3.388      0.568    224.525 

> coef(glmfit, s=glmfit$lambda.min)
59 x 1 sparse Matrix of class "dgCMatrix"
                                          1
(Intercept)                   -1.020469e+03
n_tokens_title                 7.138452e+01
n_tokens_content               2.131744e-01
n_unique_tokens                .           
n_non_stop_words               .           
n_non_stop_unique_tokens       .           
num_hrefs                      2.471492e+01
num_self_hrefs                -4.379362e+01
num_imgs                       9.106714e+00
num_videos                     3.835288e+00
average_token_length          -2.218920e+02
num_keywords                   5.555901e+01
data_channel_is_lifestyle     -2.802095e+02
data_channel_is_entertainment -8.463057e+02
data_channel_is_bus           -1.060347e+02
data_channel_is_socmed         .           
data_channel_is_tech           .           
data_channel_is_world         -1.057383e+02
kw_min_min                     2.126980e+00
kw_max_min                     2.320646e-02
kw_avg_min                     .           
kw_min_max                    -2.062232e-03
kw_max_max                    -4.173264e-04
kw_avg_max                     .           
kw_min_avg                    -2.429931e-01
kw_max_avg                    -1.435964e-01
kw_avg_avg                     1.344089e+00
self_reference_min_shares      2.098841e-02
self_reference_max_shares      3.318576e-03
self_reference_avg_sharess     .           
weekday_is_monday              4.971638e+02
weekday_is_tuesday            -6.837821e+01
weekday_is_wednesday           2.858827e+01
weekday_is_thursday            .           
weekday_is_friday             -1.962620e+01
weekday_is_saturday            4.981667e+02
weekday_is_sunday              .           
is_weekend                     1.042575e+02
LDA_00                         1.211785e+02
LDA_01                         .           
LDA_02                        -7.142931e+02
LDA_03                         8.249926e+02
LDA_04                         .           
global_subjectivity            2.209393e+03
global_sentiment_polarity      .           
global_rate_positive_words    -5.516097e+03
global_rate_negative_words     .           
rate_positive_words            .           
rate_negative_words            .           
avg_positive_polarity         -7.201715e+02
min_positive_polarity         -1.324023e+03
max_positive_polarity          .           
avg_negative_polarity         -1.340245e+03
min_negative_polarity          .           
max_negative_polarity         -6.537973e+02
title_subjectivity             .           
title_sentiment_polarity       6.161529e+01
abs_title_subjectivity         5.350562e+02
abs_title_sentiment_polarity   5.216084e+02

> coef(cvfit)[which(coef(cvfit)!=0)]
                  (Intercept)                n_tokens_title              n_tokens_content 
                -1.068705e+03                  7.338796e+01                  2.246641e-01 
              n_unique_tokens                     num_hrefs                num_self_hrefs 
                 3.767190e-01                  2.477769e+01                 -4.510284e+01 
                     num_imgs                    num_videos          average_token_length 
                 9.177435e+00                  4.384245e+00                 -2.231391e+02 
                 num_keywords     data_channel_is_lifestyle data_channel_is_entertainment 
                 5.639958e+01                 -3.061264e+02                 -8.549330e+02 
          data_channel_is_bus         data_channel_is_world                    kw_min_min 
                -1.375412e+02                 -1.168014e+02                  2.123494e+00 
                   kw_max_min                    kw_min_max                    kw_max_max 
                 2.442929e-02                 -2.107366e-03                 -4.495781e-04 
                   kw_min_avg                    kw_max_avg                    kw_avg_avg 
                -2.536708e-01                 -1.485157e-01                  1.369138e+00 
    self_reference_min_shares     self_reference_max_shares             weekday_is_monday 
                 2.104662e-02                  3.374723e-03                  5.085032e+02 
           weekday_is_tuesday          weekday_is_wednesday             weekday_is_friday 
                -6.934615e+01                  3.874155e+01                 -2.178861e+01 
          weekday_is_saturday                    is_weekend                        LDA_00 
                 5.048332e+02                  1.122331e+02                  1.643045e+02 
                       LDA_02                        LDA_03           global_subjectivity 
                -7.036686e+02                  8.096602e+02                  2.273834e+03 
   global_rate_positive_words         avg_positive_polarity         min_positive_polarity 
                -5.851422e+03                 -7.756754e+02                 -1.345198e+03 
        avg_negative_polarity         max_negative_polarity      title_sentiment_polarity 
                -1.331373e+03                 -6.922635e+02                  7.144942e+01 
       abs_title_subjectivity  abs_title_sentiment_polarity 
                 5.576846e+02                  5.330415e+02 

> library(Metrics)
> Xtest.bm <- as.big.matrix(as.matrix(d_test[,-59]))
> rmse(d_test[,59], predict(glmfit, newx=as.matrix(d_test[,-59]), s=glmfit$lambda.min))
[1] 8021.108
> rmse(d_test[,59], predict(cvfit, Xtest.bm))
[1] 8021.885
> plot(glmfit)
> plot(cvfit, type='all')

f:id:TJO:20170227143346p:plain
f:id:TJO:20170227143422p:plain


うーん。。。ワーカが立ち上がるまでにアホほど時間がかかるのも納得がいかないし、挙げ句の果てにRMSEでも{glmnet}に僅かに負けているというのが気になります。。。単に自分のMac Book Proの設定の問題かもしれません。


ちなみにご覧の通り、偏回帰係数は{biglasso}でも{glmnet}でもほとんど変わらない感じです。


特徴量選択ルールごとの違いを見てみる


一応、helpを見ると線形回帰にはデフォルトの"SSR"以外にも3種類の特徴量選択ルール(screening rule)が選べると書いてあるので、それぞれ試してみました。

> t <- proc.time()
> cvfit <- cv.biglasso(X.bm, y, family='gaussian', seed=71, ncores=3, screen = 'SEDPP')
> proc.time() - t
   ユーザ   システム       経過  
     2.985      0.605    201.546 
> rmse(d_test[,59], predict(cvfit, Xtest.bm))
[1] 8021.885

> t <- proc.time()
> cvfit <- cv.biglasso(X.bm, y, family='gaussian', seed=71, ncores=3, screen = 'SSR-BEDPP')
> proc.time() - t
   ユーザ   システム       経過  
     3.011      0.604    195.615 
> rmse(d_test[,59], predict(cvfit, Xtest.bm))
[1] 8021.885

> t <- proc.time()
> cvfit <- cv.biglasso(X.bm, y, family='gaussian', seed=71, ncores=3, screen = 'SSR-Dome')
> proc.time() - t
   ユーザ   システム       経過  
     2.774      0.597    204.556 
> rmse(d_test[,59], predict(cvfit, Xtest.bm))
[1] 8021.885

ほとんど変わらないようですorz またデータセットの内容が変わったり、データセット容量が増えたりしたら変わってくるのかもしれません。


補足


これは未だによく分かっていないんですが、helpにも載っているcolonデータを使った例だと{biglasso}のワーカは遥かに早く立ち上がってサクッと回るんですよね。何故今回使ったonline newsデータセットではアホほど時間がかかったのかは現在調査中です。。。