読者です 読者をやめる 読者になる 読者になる

琥珀色呑んだくれ備忘録

メモとか備忘録とか

Lassoの非ゼロ変数を代替する候補を可視化する

lasso feature selection interpretability R

最近、モデルの解釈性(interpretability)について色々と研究が出ており、興味もあってぼんやりと追いかけている。

“Finding Alternate Features in Lassos”という論文が提案されている: 
Satoshi Hara and Takanori Maehara (2016): "Finding Alternate Features in Lasso", in NIPS 2016 workshop on Interpretable Machine Learning for Complex Systems, December 10th, 2016. (arXiv, github
著者の日本語ページにもコンパクトな解説が掲載されている。

Lassoには変数間の相関が高いとき、それらの変数の中から一つだけを選択したがる性質がある(「見落としの問題」と呼ぶそうだ)。この問題に対して、この論文では、LASSOによって選択された変数がもし使われなかったとしたら、かわりにどの変数が選択されるかと考える。著者らによると、代替となる変数を列挙提示することで、納得感や結果の解釈性が得られるのではないか、と期待している。顧客に選択肢を残しておくというのは良いんじゃないだろうか。

具体的な手続きはシンプルで、Lassoで選択された係数\beta^*の推定量が非ゼロの変数X_iと推定量ゼロの変数X_jの組を考え、\beta^*のもとで全ての組(i,j)について\beta_i=0\beta_j\not=0とおいて目的関数を解きなおす。\beta_j\not=0と再推定されたら、「X_jX_iの代替となる変数である」と提示され、やっぱり\beta_j=0なら「X_jX_iの置き換えにはならない」とする。
また、代替となる変数同士の優劣は、目的関数の変化量を指標スコアとして相対的に順位づけされる。

で、ざっと論文を読んだ後、著者のPythonコードを見たらRでも簡単にできそうだったので、シンプルなOLS + soft-thresholdな罰則項という最小限だけをRで実装した(github)。

www.slideshare.net
ローカルに持ってきてDESCRIPTIONを用意してビルドすればパッケージになる(はず)。二分グラフを縦にプロットする方法で地味に手こずったが、この論文の可視化以外にもどこかで役に立つだろう。

面白いなと思った点は、目的関数をタスクに合わせてフレキシブルに設計できるところ。実際に、論文では損失関数をloglossに差し替えて二値分類タスクに適用・実験している。正則化項のほうを弄れば、例えばadaptive Lassoなどへの拡張もできるのではないかと思う。

また、各変数の相関行列と代替性の指標スコアとを比較すると、必ずしも大きさの順番がが一致しないのも興味深い。相関の高い変数が上位に来るのはわかりやすい一方、罰則の強さ(\lambda)の選び方によっては、割と相関の高くない変数でも代替候補として挙げられてくる。これはもしかしたら、前処理としてマルチコを除く作業をこの手法によってモデル構築後の二次工程に回せるかもしれないし、データによっては発見的な議論につながるかもしれない。
\lambdaの選び方について、スライドの例ではcross-validationで求めた値を使った。しかし、変数が増えていくと見なければいけない情報が多すぎるので、強めの罰則をかけてエッセンシャルな変数だけに絞って分析するほうがシンプルかな、とは思った。この辺りは好みもあるかもしれないが、元の論文中では推奨値に関する議論は特に言及がなかったと思う。

代替候補の変数は必ずオリジナルの変数よりも目的関数の評価が下がることもあり、おそらく予測はオリジナルの\beta^*のモデルでやったほうがよい。そういう理由でprediction()は実装しなかった。この方法で予測したいとか2値分類について同様の分析がしたい方は本家のPython版を利用されたい。

推定値がズレる縮小推定の問題や、サンプルサイズよりも変数の数が大きいときに見過ごすといった、Lassoの癖というか性質に対処するアプローチの一例については、以前紹介したことがあるのでそちらも参照されたい。

www.slideshare.net

mxnetパッケージをUbuntuで使うためのセットアップ

R deep learning mxnet ubuntu install_deps

Ubuntu14.04でRからmxnet使いたくなってセットアップ。ガイド見ながらやっていけば終わり...と思ったらinstall_depsでいくつか引っかかった。Windows/Macばりに甘やかされたい。

次からコピペでやりたいので自分用にメモしておく。
http://mxnet.readthedocs.io/en/latest/how_to/build.html#r-package-installationmxnet.readthedocs.io

# install mxnet
[~]$ sudo apt-get install build-essential git libatlas-base-dev libopencv-dev
[~]$ git clone --recursive https://github.com/dmlc/mxnet
[~]$ cd mxnet
[~/mxnet]$ make

今回、makeは特に失敗せず完了。makeがおわったら以下を実行する。インストールガイドどおりにコマンドラインからやっても良い。

require(devtools); 
require(methods);
require(httr)

setwd("mxnet/R-package/")
set_config(use_proxy(url="プロキシ使ってれば指定",port=8080))
options(repos=c(CRAN='https://cran.rstudio.com'));
install_deps(dependencies = TRUE)

install_depsが動作するディレクトリに注意。今回はプロキシ指定とかが(なぜか)上手くいかなかったので、コマンドラインじゃなくスクリプトであれこれ試行錯誤した。

いろいろ入れてる途中でimagerでfftw3.hがないよと言われて失敗したので追加。

[~]$ sudo apt-get install fftw3 fft3-dev
install_depsが沈黙したら下記を実行。
>|sh|
[~]$ cd mxnet
[~/mxnet]$ make rpkg

あとはインストール。

[~]$ R CMD INSTALL mxnet_0.7.tar.gz

以上をシェルスクリプトにまとめれば終わり。

random forestの感度分析+

R random forest mlR

 教科書を読むと、random forestでモデルを学習させて、ある良い予測精度のものが得られたら、まずは各変数の重要度をみて、つぎにpartial dependency plot(PDP)を眺めなさい、と書いてある。ある変数の動きに注目して予測値がどう変化するか見てみよう、という分析だ。

 これを発展(?)させたforestFloorは、決定木の各ノードでの予測値の変化量に着目して、random forestにおける変数が動いたとき予測値がどちら向きにどのくらい動くか定量化しようというもの。 

www.slideshare.net

スライド中のコードは、ほとんどマニュアルの丸写しなので割愛する。

 forestFloorでは、ある決定木でのある事例の予測値(もしくは、あるクラスへの所属確率)= ブートストラップ平均 + ノード1での変化量 + ...  + 終端ノードの変化量、と解釈する。この解釈の下で、ブートストラップによって当該の観察事例が出てこない木や、観察事例が終端ノードに到達するまでに注目する変数を使わなかった、というケースは除外して、注目する変数による予測値の変化量をすべての木について足し合わせることで、ある観察事例に対するある変数の貢献度(feature contribution)を計算する。

 このアイデアが面白いのは、回帰(→予測値)であろうが分類(→あるクラスへの所属確率)であろうが、弱学習器に決定木を使っているアンサンブル手法であれば貢献度が計算できる点だと思う。 弱点は、論文を眺めた限りでは、他の変数の影響を排除できているのかよく判らない点と、「それで結局予測はどうなるんだ」と訊かれそうな点と、原理上アンサンブルツリー縛りがある点。

 とくに後者2つで、一般性(や有用性)という観点からPDP&ICEboxに譲ってしまうと思われる。Hastieたちの教科書では、partial dependenceは任意の教師あり学習に対して、ある変数についてのmarginal averageが計算できると説明されている。これに加えて、ICEBoxを使うと予測値がの動きを分解(disaggregate)でき、かつ、予測値が大きく変化する変数値のポイントを検出できる。一応、forestFloorでも似たようなことはできるし、論文でもICEboxについて言及してはいるのだが、若干歯切れの悪い流れであった。。

 ちなみに、これらはmlRにて任意の教師あり学習に適用可能なwrapperとして実装されていて、下記のリンク先ではSVMをつかった謎のデータirisの3クラス分類での例が示されている。

Partial Prediction Plots - mlr tutorial

 

(追記)

上記のリンクがいろいろと切れてると思ったら、mlr 2.8→2.9で関数名の変更があったようだ。

Partial Dependence Plots - mlr tutorial

LassoとOracle Property

Oracle Property」。巷にこれほどかっこいいフレーズはそうそうないと思うので紹介を試みる。

さて、Lassoの罰則項はどういう作用を持っているか? L1 normの縮小推定を解説する図はいろんな教科書にも載っているのだが、自分は最初見たときに実感を持てなかった。Rで実際に条件を変えながらOLSと比較することで、一助になれば幸いである。

www.slideshare.net


今回使用したコード。 foreachのブロックは後半で結構時間がかかるので、イテレーションを減らすか並列化されたい。

# comparison among lm, lasso, adaptive lasso, and rigorous lasso

require(dplyr)
require(foreach)

#install.packages("glmnet")
require(glmnet)
#install.packages("msgps")
require(msgps)
#install.packages("hdm")
require(hdm)


# functions ---------------------------------------------------------------

# generate data according to model: y ~ 0 + B_i*x_i + N(0,1) 
Data_with <- function(n,p,s, silent=TRUE){
  # data with sample size(n), number of variables(p) and non-zero vector(s)
  if(! silent){
    cat(paste0(n, " obs. of ", p, " variables.\n",
              "Last ", length(s)," variables are non-zero."))
  }
  beta <- c(rep(0,p-length(s)),s)
  X <-  matrix(rnorm(n*p), ncol=p)
  
  res <-list(
    X = X,
    y = 0 + X%*%beta + rnorm(n),
    beta = beta,
    non.zero.coeff = which(beta>1e-8),
    non.zero.value = s 
  )
  if(! silent){
    str(res)
  }
  invisible(res)
}


# effect of shrinkage on estimates (lm vs lasso) ----------------------------------

test_data <- Data_with(n=10^3, p=10, s= c(1,10,100))
str(test_data,max.level = 2)
attach(test_data)

system.time(  lm.reg <- lm(y~X)  )
summary(lm.reg)

system.time(  lasso_cv  <- cv.glmnet(x = X, y=y, alpha = 1)  )
coef(lasso_cv, s="lambda.1se")

test_data <- Data_with(n=10^3, p=10, s= c(100,100,100))
str(test_data,max.level = 2)
attach(test_data)

system.time(  lm.reg <- lm(y~X)  )
summary(lm.reg)

system.time(  lasso_cv  <- cv.glmnet(x = X, y=y, alpha = 1)  )
par(mfrow=c(1,2))
plot(lasso_cv)
abline(v=log(lasso_cv$lambda.1se), lty=2, col="blue")

# plot solution path (Lasso) ------------------------------------------------------
lasso.reg <- glmnet(x = X, y=y, alpha = 1) 
plot(lasso.reg, xvar="lambda", label=TRUE)
abline(v=log(lasso_cv$lambda.1se), lty=2, col="blue")
abline(h= c(1,2,3), lty=2, col= c(1,2,3))
par(mfrow=c(1,1))

fit.coeff <- coef(lasso_cv, s="lambda.1se") %>% summary %>% transmute(beta=i-1, Var=x)
coef(lasso_cv, s="lambda.1se")
fit.coeff %>% str

# using lasso with oracle property (adaptive lasso) -------------------------

system.time(  alasso.reg <- msgps(X= X, y=c(y), penalty="alasso",gamma=1, lambda=0)  )

par(mfrow=c(2,2))
plot(alasso.reg, stand.coef=FALSE, criterion="cp",  main="Mallows' Cp")
abline(h= c(1,2,3), lty=2, col= c(1,2,3))

plot(alasso.reg, stand.coef=FALSE, criterion="aicc",main="bias-corrected AIC (AICc)")
abline(h= c(1,2,3), lty=2, col= c(1,2,3))

plot(alasso.reg, stand.coef=FALSE, criterion="gcv", main="generalized cross validation (GCV) ")
abline(h= c(1,2,3), lty=2, col= c(1,2,3))

plot(alasso.reg, stand.coef=FALSE, criterion="bic", main="BIC")
abline(h= c(1,2,3), lty=2, col= c(1,2,3))
par(mfrow=c(1,1))
summary(alasso.reg)
ss=300

# measure: estimate error vs sample size -------------
sample_sizes <- c(5:20*10, 400, 600, 800, 1000,2000,3000, 4000)
itr <- 1:10
conditions <- expand.grid(ss=sample_sizes, itr=itr) %>% arrange(ss)

sse.sampsi <- foreach(Ci=1:NROW(conditions), .combine = rbind) %do% {
  
  cond <- conditions[Ci,] 
  this_data <- Data_with(n=cond$ss, p=10, s= c(10,10,10), silent=TRUE)
  beta <- this_data$beta
  
  fit.lm <- lm(this_data$y~this_data$X)
  sse.lm <- sum( (beta - coef(fit.lm)[-1])^2 )
  
  fit.lasso_cv   <- cv.glmnet(x = this_data$X, y=this_data$y, alpha = 1)
  coeff.lasso_cv <- coef(fit.lasso_cv, s="lambda.1se") %>% 
    summary %>%
    transmute(beta=i-1, Var=x)
  sse.lasso_cv <- sum( (beta[coeff.lasso_cv$beta[-1]] - coeff.lasso_cv$Var[-1])^2 )
  
  fit.alasso   <- msgps(X= this_data$X, y=c(this_data$y), 
                        penalty="alasso", gamma=1, lambda=0)
  coeff.alasso <- coef(fit.alasso) %>% 
    data.frame %>% 
    select(BIC) %>%
    unlist
  sse.alasso   <- sum( (this_data$beta - coeff.alasso[-1])^2 )
  
  data.frame(sample.size=cond$ss, lm=sse.lm, lasso=sse.lasso_cv, alasso=sse.alasso)
}

sse.summ.sampsi <- sse.sampsi %>% 
  group_by(sample.size) %>%
  summarize_each(funs(mean))
sse.summ.sampsi

with(sse.summ.sampsi,{
  plot (lm   ~sample.size, type="l", lwd=2, col="green", ylim=c(0,0.35), ylab="mean SSE", log="x")
  lines(lasso~sample.size, type="l", lwd=2, col="red")
  legend("topright", legend = c("lm","lasso"), col = c("green","red"), lty = 1, lwd=2)
})

with(sse.summ.sampsi,{
  plot (lm    ~sample.size, type="l", lwd=2, col="green", ylim=c(0,0.35), ylab="mean SSE", log="x")
  lines(lasso ~sample.size, type="l", lwd=2, col="red")
  lines(alasso~sample.size, type="l", lwd=2, col="blue")
  legend("topright", legend = c("lm","lasso","a-lasso"), col = c("green","red","blue"), lty = 1, lwd=2)
})

# example(rlasso) ---------------------------------------------------------
# use Lasso with rigorous penalty
system.time( 
  nonpost_lasso.reg <- rlasso(
    formula, data, 
    post      = FALSE, 	# FALSEで、罰則だけrigorousなLasso
    intercept = TRUE,	  # TRUEで切片項(β0)の罰則付き推定をしない
    penalty = list(homoscedastic = FALSE, 	     # FALSEで誤差の等分散性を仮定しない
                   X.dependent.lambda = FALSE,   # TRUEでデータ依存的なλの選択をする
                   lambda.start = NULL, 
                   c = 1.1, 
                   gamma = 0.1/log(n)),
    control = list(numIter =15,
                   tol = 10^-5,
                   threshold = NULL),
    ...)
  
)

# use Post-Lasso procedure with rigorous penalty
system.time( 
  post_lasso.reg <-  rlasso(
    formula, data, 
    post      = TRUE, 	# TRUEでPost-Lasso+rigorousな罰則
    intercept = TRUE,	  # TRUEで切片項(β0)の罰則付き推定をしない
    penalty = list(homoscedastic = FALSE, 	     # FALSEで誤差の等分散性を仮定しない
                   X.dependent.lambda = FALSE,   # TRUEでデータ依存的なλの選択をする
                   lambda.start = NULL, 
                   c = 1.1, 
                   gamma = 0.1/log(n)),
    control = list(numIter =15,
                   tol = 10^-5,
                   threshold = NULL),
    ...)
  
)
# comparison results (set all=FALSE to summarize)
summary(nonpost_lasso.reg, all=FALSE)
summary(post_lasso.reg, all=FALSE)
coef(post_lasso.reg)


# measure: estimate error vs dimension size -------------
dim_sizes <- c(1:10*10, 1:5*200, 2:5*1000, 1e4, 2e4)
itr <- 1:30
conditions <- expand.grid(ds=dim_sizes, itr=itr) %>% arrange(ds)

sse.dim_sizes <- foreach(Ci=1:NROW(conditions), .combine = rbind) %do% {
  
  cond <- conditions[Ci,] 
  this_data <- Data_with(n=1000, p=cond$ds, s= c(10,10,10), silent=TRUE)
  beta <- this_data$beta
  
  fit.lm <- lm(this_data$y~this_data$X)
  sse.lm <- sum( (beta - coef(fit.lm)[-1])^2 )
  
  fit.lasso_cv   <- cv.glmnet(x = this_data$X, y=this_data$y, alpha = 1)
  coeff.lasso_cv <- coef(fit.lasso_cv, s="lambda.1se") %>%
    summary %>% 
    transmute(beta=i-1, Var=x)
  sse.lasso_cv   <- sum( (beta[coeff.lasso_cv$beta[-1]] - coeff.lasso_cv$Var[-1])^2 )
  
  fit.alasso    <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1)
  coeff.alasso  <- coef(fit.alasso) %>% 
    data.frame %>%
    select(BIC) %>%
    unlist
  sse.alasso    <- sum( (this_data$beta - coeff.alasso[-1])^2 )
  
  plasso.reg    <- rlasso(c(this_data$y) ~ this_data$X, post=TRUE, intercept=TRUE)
  coeff.plasso  <-coef(plasso.reg)
  sse.plasso    <-(this_data$beta - coeff.plasso)^2 %>% sum
  
  data.frame(dim.size=cond$ds, 
             lm=sse.lm, lasso=sse.lasso_cv, alasso=sse.alasso, plasso=sse.plasso)
  
}

sse.summ.dim_sizes <- sse.dim_sizes %>% 
  group_by(dim.size) %>% 
  summarize_each((funs(mean)))
sse.summ.dim_sizes

with(sse.summ.dim_sizes,{
  plot( lm    ~dim.size, type="l", lwd=2, col="green", ylab="mean SSE", log="x")
  lines(lasso ~dim.size, type="l", lwd=2, col="red")
  lines(alasso~dim.size, type="l", lwd=2, col="blue")
  
  legend("topright", legend = c("lm","lasso","a-lasso"),
         col = c("green","red","blue"), lty = 1, lwd=2)
  abline(v=1000, lty=3)
})

with(sse.summ.dim_sizes,{
  plot(lasso~dim.size, type="l", lwd=2, col="red", ylim=c(0,0.1), ylab="mean SSE",log="x")
  lines(alasso~dim.size,  type="l", lwd=2, col="blue")
  
  legend("topleft", legend = c("lasso","a-lasso"), col = c("red","blue"), lty = 1, lwd=2)
  abline(v=1000, lty=3)
})


with(sse.summ.dim_sizes,{
  plot(lasso~dim.size, type="l", lwd=2, col="red", ylim=c(0,0.1), ylab="mean SSE",log="x")
  lines(alasso~dim.size,  type="l", lwd=2, col="blue")
  lines(plasso~dim.size,  type="l", lwd=2, col="orange")
  
  legend("topleft", legend = c("lasso","a-lasso","post-lasso"), 
         col = c("red","blue","orange"), lty = 1, lwd=2)
  abline(v=1000, lty=3)
})


# measure: performance vs dimension size -------------
this_data <- Data_with(n=1000, p=2e4, s= rep(10,3),silent = FALSE)

system.time(
  plasso.reg <- rlasso(c(this_data$y) ~ this_data$X, post=TRUE, intercept=TRUE)
)

#coef(fit.lasso_cv)
system.time(
  fit.alasso <- msgps(X= this_data$X, y=c(this_data$y), penalty="alasso", gamma=1)
)


反省点その1、これほどかっこいいフレーズなのに当日聴衆には全く刺さらなかった(というか完全に滑った)。プレゼンスキルの改善がである。その2、圧縮センシングのお得感とか、スパース推定後の解釈性とか、もうちょっと使い道に触れた導入にすべきだった。最後の"rigorous lasso"では、計算時間は少なからず増えているわけで、実際の現場感としてここまでの精度が必要か?とか、いずれ費用対効果やロバスト化の議論はしてみたい。

random forestを予測以外の目的で使う

R random forest

数年ほど前には最強と言われて一世を風靡していたrandom forestだが、予測以外にも使い道が提案されている。Rのパッケージから紹介したい。

 予測全体の把握と仮説ルールの抽出

決定木分析が便利な理由の一つは「どういうルールでその予測が成り立っているのか」を極めて簡単に可視化出来る点。inTreesパッケージは予測ルールを集計し、適当に枝狩りして全体を要約することで、アンサンブルモデルにおいても決定木と同じような情報を可視化してくれる。

www.slideshare.net

このパッケージではもう一つ、すべての木から取り出した1つずつの枝をトランザクションとみなしてアソシエーション分析する機能を提供する。メジャーな関連ルールを取り出す方に力点が置かれているが、変数同士の稀な関係が意外なアウトカムを導いてる方が、(研究上は)見つかったときに面白いので、欲を言えばLift値も評価してくれると良かった。

ノンパラメトリックな欠損値の補完

randomForestパッケージには、rfImputeという欠損値を補完する関数が用意されている。ところが、rfImputeは同じルールにより同じ予測をされる別の観察データを使って補完するため、目的変数が欠損しているデータには適用できない。

既知のデータにおいて欠測部分を推定するにはこれで構わないのだが、予測をしたいデータの前処理に使おうとするとこれでは困る。

missForestでは、補完対象の変数を目的変数とみなすことで、random forestにより欠損値を直接予測する。

www.slideshare.net

もともと予測したい変数も補完プロセスに含めることができるので半教師学習と言える。その場合、イテレーションの回数と過学習の関係が気になるが、元の論文にはその言及はなかったように思う。

 

“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%をネストしたい

R

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

評価関数は例えば、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

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