琥珀色呑んだくれ備忘録

メモとか備忘録とか

“rigorous” lassoを試してみた

High-Dimensional Metrics in R
CRAN - Package hdm

  • やりたかったこと= oracle property保証のlassoで変数選択(n<p)
  • 得たもの= 速い、一致性はそれっぽい。
  • やりたいこと= ほかのパッケージとの比較、証明部分読む。

以下、試したかった部分だけ:

# install.packages("hdm")
require(hdm)
set.seed(1)

### modified from example(rlasso) ----------------

n <- 100  # sample size
p <- 1500 # number of variables
s <- 3    # nubmer of non-zero variables
beta = c(rep(3,s), rep(0,p-s)) # s -> non-zero/zero coeff.vec

X <- matrix(rnorm(n*p), ncol=p)
y <- 1 + X%*%beta + rnorm(n)

# use Lasso,not-Post-Lasso 
system.time(
  lasso.reg <- rlasso(y~X,post=FALSE, intercept=TRUE)
)
summary(lasso.reg, all=FALSE) # use this option to summarize results

# use Post-Lasso 
system.time(
  lasso.reg <- rlasso(y~X, post=TRUE, intercept=TRUE)
)
summary(lasso.reg, all=FALSE) # use this option to summarize results


yhat.lasso <- predict(lasso.reg)   #in-sample prediction
Xnew <- matrix(rnorm(n*p), ncol=p)  # new X
Ynew <- Xnew%*%beta + rnorm(n)  #new Y

yhat.lasso.new = predict(lasso.reg, newdata=Xnew)  #out-of-sample prediction

オペレーションズ・リサーチ(Vol.58, No.5, 261)の解説記事にいろんなlassoとRの対応パッケージ表があるので、参考例のデータは微妙だがadaptive lassoとかとも比較してみたい。
http://www.orsj.or.jp/archive2/or58-05/or58_5_261.pdf
ちなみにこの記事、線形回帰モデルからlassoへの流れ、L1正則化が変数選択として機能する理由、変数選択のオラクルプロパティ、その先の展開まで、わずか数ページでまとまっていてオススメ。

詳しい解説
http://www.is.titech.ac.jp/~s-taiji/tmp/sparse_tutorial_2014.pdf

そのうちきちんと書き直します。

foreach+%dopar%をネストしたい

並列化できる評価関数があり、その調節パラメータを並列にグリッドサーチしたい。単純に考えると親のソケットクラスタから子クラスタを生成すると良いのかな、と思うがもっといい方法ないものか。

評価関数は例えば、randomForestのようなアンサンブル学習でも良いし、遺伝的アルゴリズムのような最適化でも良いし、LOOみたいなcross-validationを考えても良い。Rで簡単に書ける。

require(foreach)
require(doParallel)

eval_parallel <- function(){
  # 親ソケットクラスタ
  cl <- makeCluster(8)
  registerDoParallel(cl)
  on.exit(stopCluster(cl))

  # 子クラスタの格納用
  cl_i <-NULL

  # grid searchの条件組み合わせを列挙しておく
  conditions <- expand.grid(a=1:3, b=letters[2:7], c=c(1.1,1.7,1.9))

  foreach(i=1:NROW(conditions), .packages = c("foreach", "doParallel"), .combine = rbind) %dopar%{
    # 子クラスタ
    cl_i[[i]] <- makeCluster(11)
    registerDoParallel(cl_i[[i]])
    on.exit(stopCluster(cl_i[[i]]))

    # ハイパーパラメータのセット
    parameters <- conditions[i,]

    # iteration_sizeは弱学習器の数だったり、CVのfolding数だったり、GAの個体数だったり
    foreach(j=1:iteration_size, .combine = c) %dopar% {
      # do something
    }
  }
}

eval_parallel()

この例では親が8いて子が11づつ作られるので、8*(1+11)=96並列。親子のバランスはボトルネックを見て配分すればよい。ときどき親の割り当てぶんの確保を忘れる。
.combineがrbindとcになってるのは、このあとでrowMeansすれば平均精度を見てベストなパラメータセットを探せるcross-validationのイメージ。

on.exit()使ってるのは下記の記事から。
R の foreach で並列処理するとき on.exit() が便利 #rstatsj - Qiita

もうちょっとうまいやり方がありそうだけど、思いつかないので今はこれで凌いでいる。