R 3.4.0 になり, for文がデフォルトでコンパイルされることで
以前に比べて R の計算速度は早くなりました。
何よりも R の計算のボトルネックはとにかく for文でした。
（いかに for文を使わずにプログラミングできるかが
結構重要なスキルだった時代も、、、・）
そんな高速化している R ですが、
iteration × Cross validation が当たり前になるなど
それ以上の速度で我々が扱う計算は
タスクが大きくなってます。
ゆえに相変わらず重要なのが
「並列化」。
特に、最近のPCでは10個近い
並列計算が簡単になってきていて
これを使わない手はありません。
というわけで
以前の個人HPでも紹介していた
Rで “簡単に” 並列計算をする方法を
改めて紹介します。
簡単に内容をまとめると、、、
・foreach と doParallel というパッケージを使うと
並列計算が超簡単
・何が簡単かというと、
「各パターンごとに計算」というのが
シンプルに実行できる
・i.e. Monte Carlo シミュレーションには最適
という感じ。
というわけで実装してみました。
プログラム内容は
・coordinate descent algorithm (glmnetパッケージ)
・elastic net
・10-fold cross-validation (alpha と lambda の grid 検索)
です。
ちなみに glmnet パッケージも進化しており、
最近ではデフォルトで iterators パッケージと
foreach パッケージを読み込むようになりました。
【準備：CRAN から package をゲット】
1. glmnet：https://cran.r-project.org/web/packages/glmnet/index.html
（elastic net の計算で使用. 並列計算とはあまり関係なし）
2. iterators：http://cran.r-project.org/web/packages/iterators/index.html
（下の foreach のために必要らしい. glmnet を使うなら不要）
3. foreach：http://cran.r-project.org/web/packages/foreach/index.html
（glmnet を使うなら不要）
4. doParallel：http://cran.r-project.org/web/packages/doParallel/index.html
【プログラム】
library(glmnet)
library(MASS)
library(iterators)
library(foreach)
library(doParallel)
registerDoParallel(detectCores())
### 人工データ作成
n <- 50
p <- 100
sigma <- 5
sigma2 <- sigma**2
rho <- 0.5
t_beta <- c(rep(3, ceiling(p/2)), numeric(p-ceiling(p/2)))
Sigma <- diag(p)
for(i in 1:p){
for(j in 1:p){
Sigma[i,j] <- rho**abs(i-j)
}
}
x <- mvrnorm(n, numeric(p), Sigma)
H <- diag(n) - matrix(rep(1/n, n**2), nrow=n)
x <- H %*% x
x <- x %*% diag(sqrt(n)/sqrt(colSums(x**2)))
y <- x %*% t_beta + rnorm(n, 0, sigma)
y <- y-mean(y)
### 候補の lambda と alpha の数の決定（今回は 100×100 の grid search）
lg.lam <- 100
lg.al <- 100
### k-fold cross-validation の k
K <- 10
### alpha fix での k-fold cross-validation
kcv.en <- function(x, y, alpha, lg.lam, K){
n <- nrow(x)
p <- ncol(x)
y <- drop(y)
lambda.max <- max(abs(t(x) %*% y))/(n*alpha)
lambda.min <- 0.01
lambdas <- (exp(((log(lambda.max)-log(lambda.min))
* ((lg.lam-1):0)/(lg.lam-1)) + log(lambda.min)))
PSE <- matrix(numeric(lg.lam*K), nrow=K)
n1 <- n
N <- matrix(numeric(2*n), nrow=n)
Ind <- sample(1:n, n, rep=FALSE)
N[,1] <- Ind
for(k in 1:K){
n2 <- ceiling(n1/(K-k+1))
N[(1:n2)+(n-n1),2] <- k
n1 <- n1 - n2
}
for(k in 1:K){
X.est <- x[N[,2] != k, ]
y.est <- y[N[,2] != k]
X.evl <- x[N[,2] == k, ]
y.evl <- y[N[,2] == k]
GLMNET <- glmnet(X.est, y.est, family=\"gaussian\", alpha=alpha, lambda= lambdas)
Beta <- t(matrix(as.numeric(GLMNET$beta), nrow=p))
y.prd <- X.evl %*% t(Beta)
PSE[k,] <- colSums((y.prd - y.evl)**2)
}
CV <- colSums(PSE)/K
### CV の値と lambda の値を返す関数
return(list(cv=CV, lambda=lambdas))
}
alphas <- (1:lg.al)/lg.al
CV <- matrix(numeric(lg.lam*lg.al), ncol=lg.lam)
LAMBDA <- matrix(numeric(lg.lam*lg.al), ncol=lg.lam)
### 並列計算なし
TIME1 <- system.time(for(i in 1:lg.al){
alpha <- alphas[i]
res <- kcv.en(x, y, alpha, lg.lam, K)
CV[i,] <- res$cv
LAMBDA[i,] <- res$lambda
})
image(CV)
### 並列計算あり
f.kcv <- function()
{
foreach(i=1:lg.al, .export=ls(envir=parent.frame())) %dopar% {
alpha <- alphas[i]
kcv.en(x, y, alpha, lg.lam, K)
}
}
f.para <- function(lg.lam, lg.al){
res <- f.kcv()
CV.p <- matrix(numeric(lg.lam*lg.al), ncol=lg.lam)
LAMBDA.p <- matrix(numeric(lg.lam*lg.al), ncol=lg.lam)
for(i in 1:lg.al){
CV.p[i,] <- res[[i]]$cv
LAMBDA.p[i,] <- res[[i]]$lambda
}
return(list(CV=CV.p, LAMBDA=LAMBDA.p))
}
TIME2 <- system.time(RES <- f.para(lg.lam, lg.al))
以上のプログラムで
並列計算ありとなしで時間を計測したところ、
・並列なし：4.306
・並列あり：1.234
とだいたい3〜4倍速いという結果に！！！
（ちなみにPC環境はコア数8）
が、いくつかの場面でこの方法を使用したところ
以下の注意点に気づきました
1. 必ずしも並列計算の方が速いわけではない
：ジョブを分割するためか、
実際の計算以外の時間もかかるらしく、
簡単なプログラムでは並列なしの方が速いことも。
2. メモリオーバーに注意！
：大きなデータを扱う場合、
ジョブを分割した分だけメモリが必要に。
その結果、メモリオーバーを起こしたのか
エラーが返ってくることが何度かありました。
というわけで今回は私でもできる
Rでの簡単並列計算の方法を紹介致しました。
また、並列計算と同じく
「複数マシン」を使うものとして
分散処理がありますが
最近では Hadoop以外にも選択肢があるようで、
Rとの連携という意味では Spark が
中々便利そうです。
（Spark についてもは今後勉強していくので
また機会があれば共有します。）
みなさんもRで並列計算をガンガン使ってみてくださいな。
参考文献：
http://d.hatena.ne.jp/teramonagi/20140920/1411195147
