今回は小ネタ。ボサーっと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')
うーん。。。ワーカが立ち上がるまでにアホほど時間がかかるのも納得がいかないし、挙げ句の果てに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データセットではアホほど時間がかかったのかは現在調査中です。。。