(2023S) 1152-4001 Statistical consulting

*email address: summet73y.da@gmail.com

Faster computations in R

• Topic 4에서는 R을 통한 계산을 보다 빨리 수행할 수 있는 다음의 방법들에 대해서 다룰 것이다.
  1. 빠른 R 계산을 위한 여러가지 팁
  2. 병렬 계산 (Parallel computing)
  3. Rcpp 패키지 사용법.
• R에서는 보다 빠른 계산 및 실행을 위하여 C나 C++을 기반으로 한 패키지를 활용하는 경우가 많다.
  • Windows 환경의 경우, C나 C++을 기반으로 한 패키지를 보다 안정적으로 구동하기 위하여 R버전에 맞게 Rtools를 추가적으로 설치해야한다.

Section 1. 빠른 R 계산을 위한 팁

1.1 Bottleneck

• R코드를 빠르게 실행하기 위해서는 R코드 상에서 실행하는데 가장 큰 시간이 소모되는 Bottleneck을 효율적으로 바꾸는 것이 가장 중요.
• 즉, 다음과 같은 과정을 통해 R 코드를 효율적으로 바꿀 수 있다.
  1. R코드 상에서 Bottleneck 인지
  2. Bottleneck 부분을 보다 빠르게 실행시킬 수 있도록 효율적으로 코드를 다시 작성.
• 이 때, 1번 과정에서 Bottleneck을 인지하는 것은 Code profiling이라고도 한다.
  • Code profiling은, 작성된 코드의 구조에서 각 부분에서의 메모리와 시간이 어느정도 소모가 되는지를 분석하는 작업을 뜻한다.
• 실제로 작성된 코드에서 Bottleneck을 인지하는 것은
  • 일반적으로는 어느 부분이 가장 큰 시간을 소모하는지 코드 작성자가 알지만,
  • 다른 이가 작성한 코드를 사용하거나, 본인이 작성한 코드라고 하더라도 코드가 복잡한 경우에는 어느 부분이 Bottleneck인지 모르는 경우가 있을 수 있다.
• Bottleneck 부분이 어느 곳인지 모르는 경우는, 이를 확인하기 위해 실제로 R코드를 실행하면서 각 부분에서 걸리는 시간을 확인하면서 진행한다.
  1. system.time() 함수 사용
  2. tictoc 패키지 활용
  3. Rprof() 함수 사용
  4. profvis 패키지 사용

system.time() 함수

###########################################
# Method 1 - system.time()
###########################################

### library
library(ggplot2movies)
library(profvis)

### system.time
system.time({ data(movies, package = "ggplot2movies")  })

 사용자  시스템 elapsed 
      0       0       0 

system.time({ plot(movies$year, movies$rating)  })

 사용자  시스템 elapsed 
   0.00    0.06    0.42 

system.time({ model <- loess(rating ~ year, data = movies)  })

 사용자  시스템 elapsed 
   5.06    0.02   10.35 

system.time({ j <- order(movies$year) })

 사용자  시스템 elapsed 
      0       0       0 

system.time({ lines(movies$year[j], model$fitted[j]) })

 사용자  시스템 elapsed 
   0.00    0.00    0.01 

• 결과에서
  • 사용자 : CPU에서 코드를 실행하는데 걸리는 시간,
  • 시스템 : 시스템 프로세스에서 걸리는 시간. 예를 들어 파일을 열거나 닫는 데 걸리는 시간.
  • elapsed : 실제로 코드가 실행되는데 걸리는 시간.

tictoc 패키지

• system.time()과 수행하는 일은 거의 비슷하지만, 보다 편리하게 사용할 수 있는 tictoc패키지가 있다.

###########################################
# Method 2 - tictoc package
###########################################

### library
library(ggplot2)
library(tictoc)


### system.time
tic("Total")
# step 1
tic("Step 1")
data(diamonds, package = "ggplot2")
toc()

Step 1: 0 sec elapsed

# step 2
tic("Step 2")
plot(price ~ carat, data = diamonds)
toc()

Step 2: 0.43 sec elapsed

# step 3
tic("Step 3")
m <- lm(price ~ carat, data = diamonds)
toc()

Step 3: 0.01 sec elapsed

# step 4
tic("Step 4")
abline(m, col = "red")

toc()

Step 4: 0 sec elapsed

toc()

Total: 0.46 sec elapsed

Rprof() 함수

• sample.interval만큼의 interval을 기준으로 해서 각 함수의 시간을 측정하는 것.
• sampling.time은 총 걸린 시간.
• by.total과 by.self의 차이점은 다음과 같다.
  • 기본적으로 R함수 안에서 다른 함수를 같이 실행하는 경우에 대해서
    • total은 하나로 합쳐서 나타내는 반면에
    • self는 함수별로 따로 나타냄.

###########################################
# Method 3 - Rprof()
###########################################

### library
library(ggplot2)

### Rprof()
Rprof("result.out", interval=0.005) # Turn on the profiler
data(diamonds, package = "ggplot2")
plot(price ~ carat, data = diamonds)
m <- lm(price ~ carat, data = diamonds)
abline(m, col = "red")

Rprof(NULL)        # Turn off the profiler
summaryRprof("result.out")   # result

$by.self
           self.time self.pct total.time total.pct
"plot.xy"      0.110    78.57      0.110     78.57
"deparse"      0.025    17.86      0.025     17.86
"sys.call"     0.005     3.57      0.005      3.57

$by.total
                          total.time total.pct self.time self.pct
"block_exec"                   0.140    100.00     0.000     0.00
"call_block"                   0.140    100.00     0.000     0.00
"eng_r"                        0.140    100.00     0.000     0.00
"eval"                         0.140    100.00     0.000     0.00
"eval_with_user_handlers"      0.140    100.00     0.000     0.00
"evaluate"                     0.140    100.00     0.000     0.00
"evaluate::evaluate"           0.140    100.00     0.000     0.00
"evaluate_call"                0.140    100.00     0.000     0.00
"handle"                       0.140    100.00     0.000     0.00
"in_dir"                       0.140    100.00     0.000     0.00
"in_input_dir"                 0.140    100.00     0.000     0.00
"knitr::knit"                  0.140    100.00     0.000     0.00
"process_file"                 0.140    100.00     0.000     0.00
"process_group"                0.140    100.00     0.000     0.00
"process_group.block"          0.140    100.00     0.000     0.00
"rmarkdown::render"            0.140    100.00     0.000     0.00
"timing_fn"                    0.140    100.00     0.000     0.00
"withCallingHandlers"          0.140    100.00     0.000     0.00
"withVisible"                  0.140    100.00     0.000     0.00
"do.call"                      0.135     96.43     0.000     0.00
"plot"                         0.135     96.43     0.000     0.00
"plot.default"                 0.135     96.43     0.000     0.00
"plot.formula"                 0.135     96.43     0.000     0.00
"plot.xy"                      0.110     78.57     0.110    78.57
"deparse"                      0.025     17.86     0.025    17.86
"deparse1"                     0.025     17.86     0.000     0.00
"paste"                        0.025     17.86     0.000     0.00
"sys.call"                     0.005      3.57     0.005     3.57
"%in%"                         0.005      3.57     0.000     0.00
".External2"                   0.005      3.57     0.000     0.00
"["                            0.005      3.57     0.000     0.00
"[.data.frame"                 0.005      3.57     0.000     0.00
"[["                           0.005      3.57     0.000     0.00
"[[.data.frame"                0.005      3.57     0.000     0.00
"lm"                           0.005      3.57     0.000     0.00
"model.frame.default"          0.005      3.57     0.000     0.00
"na.omit"                      0.005      3.57     0.000     0.00
"na.omit.data.frame"           0.005      3.57     0.000     0.00
"stats::model.frame"           0.005      3.57     0.000     0.00

$sample.interval
[1] 0.005

$sampling.time
[1] 0.14

profvis패키지

https://support.posit.co/hc/en-us/articles/218221837-Profiling-R-code-with-the-RStudio-IDE 참조.

1.2 for loop를 최대한 자제

• R은 for loop의 속도가 매우 느린편에 속하기 때문에 for loop를 자제할 수 있는 경우에는 최대한 자제하는 것이 좋음
• for loop 대신 vectorize하거나 혹은 apply 관련 함수들을 적극적으로 사용함으로써 보다 빠른 계산이 가능.

Vectorize

• (예시 1) vectorize

### library
library(tictoc)

### data
x <- seq(1:100000000)

### vectorize
tic("vectorize")
x2 <- x+1
toc()

vectorize: 0.25 sec elapsed

### for loop
tic("for loop")
x2 <- x
for(i in seq(x)){ 
  x2[i] <- x[i] + 1 
}
toc()

for loop: 2.58 sec elapsed

apply 관련 함수

• apply & sapply & lapply & tapply

### library
library(tictoc)
library(dplyr)

### data
y <- matrix(1:1000000, ncol=10000)

### apply (matrix)
tic("apply")
y2 <- apply(y,2, function(x){ mean(sqrt(x)) })
toc()

apply: 0.05 sec elapsed

### sapply (vector or list / return vector or matrix)
tic("sapply")
y2 <- sapply(1:ncol(y), function(i){ mean(sqrt(y[,i])) }, simplify = T)
toc()

sapply: 0.03 sec elapsed

tic("sapply")
y2 <- sapply(as.data.frame(y), function(x){ mean(sqrt(x)) })
toc()

sapply: 0.05 sec elapsed

### lapply (vector or list / return list)
tic("lapply")
y2 <- lapply(1:ncol(y), function(i){ mean(sqrt(y[,i])) }) %>% unlist()
toc()

lapply: 0.08 sec elapsed

### for loop
tic("for loop")
y2 <- NULL
for(i in 1:ncol(y)){
  y2 <- c(y2, mean(sqrt(y[,i])))
}
toc()

for loop: 0.17 sec elapsed

### data 2
z <- 1:1000000
f <- rep(1:10000,each=100) 
tmp.dat <- data.frame(z=z,f=f)

### tapply (by group)
tic("tapply")
y2 <- tapply(z,f,function(x){ mean(sqrt(x)) })
toc()

tapply: 0.06 sec elapsed

# ### dplyr
# tic("tapply")
# y2 <- tmp.dat %>% group_by(f) %>% summarise( x = mean(sqrt(z))) %>% pull(x)
# toc()

### for loop
tic("for loop")
y2 <- NULL
for(i in unique(f)){
  y2 <- c(y2, mean(sqrt(z[f==i])))
}
toc()

for loop: 12.92 sec elapsed

for loop 안에서 object 수정은 최대한 피하기

• R은 다른 언어와는 다르게 object를 수정을 할 때 해당 object에 바로 수정하는 것이 아닌, object의 copy를 만들고 수정하는 방식을 취하는 경우가 있음. (data.frame의 경우)
• 따라서 for loop 안에서 용량이 큰 object를 수정하게 되면, 매 iteration 마다 복사를 해야하기 때문에 시간이 오래 걸림.

### library
library(tictoc)

tic("Time 1")
x <- data.frame(matrix(runif(100*1e4), ncol = 100))
y <- data.frame(matrix(runif(100*1e4), ncol = 100))
x <- rbind(x,y)
toc()

Time 1: 0.03 sec elapsed

# for loop with data.frame
tic("Time 2 : for loop with data.frame")
x <- data.frame(matrix(runif(100*1e4), ncol = 100))
for(i in seq_along(1:100)){
  x <- rbind(x, data.frame(matrix(runif(1*1e4), ncol = 100)))
}
toc()

Time 2 : for loop with data.frame: 1.05 sec elapsed

# for loop with matrix
tic("Time 3 : for loop with matrix")
x <- matrix(runif(100*1e4), ncol = 100)
system.time(
  for(i in seq_along(1:100)){
    x <- rbind(x, matrix(runif(1*1e4), ncol = 100))
  }
)

 사용자  시스템 elapsed 
   0.11    0.06    0.45 

toc()

Time 3 : for loop with matrix: 0.53 sec elapsed

Readability vs Speed

• 현재 Topic에서는 R코드를 효율적으로 빠르게 실행시키는 것에 집중. (Speed)
• 이를 위해 apply 관련 함수들을 사용할 경우, 속도가 느린 for문을 실행하지 않으면서 빠르게 계산할 수 있으나, 이는 읽기 좋은 코드는 아닐 수 있다.
• 특히 긴 코드를 작성할수록, 코드의 가독성이 중요하기 때문에 speed를 위해 readability를 희생하는 것이 좋은 선택이 아닐 때가 있을 수 있다.
• 따라서 코드를 실행할 때 bottleneck인 부분이 아니라면, 다소 speed는 희생이 되더라도 readability를 위해 for문을 사용하는 것이 더 좋은 선택일 수 있다.

1.3 더 좋은 R package 활용

• 기본적으로 R은 속도가 느리기 때문에 advanced된 R package의 경우, 위에서 설명한 Rcpp 패키지를 활용하여 C++을 통해 컴파일하여 더 빠른 속도로 계산이 가능한 경우가 있음.

  • C++ 뿐만이 아닌, Fortran 등을 이용한 패키지 등도 있으며 순수 R만 이용한 코드보다 빠름.
  • 특히, ifelse()등 R에 기본으로 내장되어있는 함수 중에서 속도가 느린 함수는 보다 효율적인 패키지의 함수로 대체해서 사용하는 편이 좋다.

• R 패키지에 대한 매뉴얼에서 (https://cran.r-project.org) 보통 관련 설명이 있는 편.

• (예시) fastglm 패키지

  • R에 기본으로 내장되어 있는 glm보다 빠르게 실행될 수 있도록, C++기반으로 generalized linear model을 fitting할 수 있는 패키지

### library
library(fastglm)
library(tictoc)

### data generation
set.seed(1)
n <- 100000
p <- 100
x <- matrix(rnorm(n*p),n,p)
y <- sample(0:1, n, replace=T)

### glm
tic("glm : ")
glm.fit0 <- glm.fit(x, y, family=binomial())
toc()

glm : : 1.46 sec elapsed

### fastglm
tic("fastglm : ")
glm.fit1 <- fastglm(x, y, family=binomial())
toc()

fastglm : : 1.18 sec elapsed

1.5 대용량 csv파일 읽기 & 쓰기

• R에서 csv파일을 읽고 쓰기 위해 사용하는 함수로는 크게 세가지로 나눌 수 있다.
  1. R 기본 내장 함수 : read.csv() & write.csv()
  2. data.table 패키지 : read_csv() & write_csv()
  3. readr 패키지 : fread() & fwrite()
• 이 때, 이 세가지 방법에 대해서 서로 시간을 비교해보면 다음과 같다.
  • 아래 결과에서 확인할 수 있다시피, data.table을 사용했을 때 가장 빠르게 실행됨을 확인할 수 있다.

### library
library(data.table)
library(readr)

### data generation
set.seed(1)
n <- 100000
p <- 100

DAT <- data.frame(matrix(rnorm(n*p),n,p))

### write
# base
tic("(write) base : ")
write.csv(DAT, "base.csv")
toc()

(write) base : : 12.1 sec elapsed

# readr
tic("(write) readr : ")
write_csv(DAT, "readr.csv")
toc()

(write) readr : : 0.2 sec elapsed

# data.table
tic("(write) data.table : ")
fwrite(DAT, "datatable.csv")
toc()

(write) data.table : : 0.06 sec elapsed

### read
# base
tic("(read) base : ")
csv1 <- read.csv("base.csv")
toc()

(read) base : : 9.89 sec elapsed

# readr
tic("(read) readr : ")
csv2 <- read_csv("readr.csv")
toc()

(read) readr : : 0.47 sec elapsed

# data.table
tic("(read) data.table : ")
csv3 <- fread("datatable.csv")
toc()

(read) data.table : : 0.16 sec elapsed

• 또한 read.csv() 및 read_csv() 그리고 fread()함수의 차이점으로 다음을 들 수 있다.
  • 각 함수로 불러온 파일의 class를 확인해보면 다음과 같다.
    1. read.csv()로 csv파일을 불러오게되면 data frame 형태의 데이터 구조로 저장이 된다.
    2. read_csv()로 csv파일을 불러오게되면 data frame 형태의 데이터 구조로 저장이 되는데, 이 때 tibble 형태로 저장이 되며, 기본적인 data frame보다 dplyr패키지와 연계하여 보다 빠르고 효율적으로 사용할 수 있다.
    3. fread()로 csv파일을 불러오게되면 data.table 형태의 데이터 구조로 저장이 되는 것을 확인할 수 있다. data.table은 tbl_df 보다도 빠르고 효율적으로 계산할 수 있다는 장점을 가진다.

### library
library(dplyr)

### class
csv1 %>% class()

[1] "data.frame"

csv2 %>% class()

[1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame" 

csv3 %>% class()

[1] "data.table" "data.frame"

• 간단하게 data.table 형태의 데이터 구조에서의 계산 속도 향상을 확인해보면 다음과 같다.
  • 각 factor별로 mean을 계산하기 위한 코드를 다음과 같이 작성하여 데이터 구조별로 비교해볼 수 있다.
  • 다만 data.table의 경우 기본 R이나 dplyr패키지를 사용할 때와 다른 문법 구조를 가지고 있기 때문에 어느정도 진입장벽이 있는 편이다.

### data generation
set.seed(1)
n <- 10000000
x <- rnorm(n)
f <- sample(1:1000, n, replace=T)
DAT1 <- data.frame(x,f)
DAT2 <- tibble(x,f)
DAT3 <- data.table(x,f)

### base
tic("base : ")
result1 <- tapply(DAT1$x,DAT1$f,mean)
toc()

base : : 1.27 sec elapsed

### dplyr
tic("dplyr : ")
result2 <- DAT2 %>% group_by(f) %>% summarise(m=mean(x))
toc()

dplyr : : 1.03 sec elapsed

### data.table
tic("data.table : ")
result3 <- DAT3[,mean(x),by='f']
toc()

data.table : : 0.06 sec elapsed

Section 2. Parallel computing

• R코드에서 보다 빠르고 효율적으로 계산하기 위한 방법으로 (특히, for 문) 병렬 계산 (Parallel computing)을 수행할 수 있음.
• 기본적으로 R코드를 실행할 경우에 하나의 코어만을 사용해서 계산을 수행.
  • for문을 수행할 때도 하나의 코어만을 수행하는데, 각 loop가 독립적이라면 이는 비효율적.
  • for문을 수행할 때 loop 여러개를 한번에 수행하도록 병렬 계산을 수행함으로써 계산 속도를 높일 수 있음.
• 이를 위해 R에서는 parallel패키지를 사용할 수 있음.
  • 이외에도 foreach 패키지와 doParallel 패키지를 활용할 수 있음.

2.1 parallel 패키지

• R에서 기본적으로 내장되어 있는 apply 관련 함수들을 보다 빠르게 여러 개의 코어를 사용해서 수행하기 위해서 parallel패키지의 parLapply()함수 등을 활용할 수 있다.
  • Mac이나 Linux 환경에서는 mclapply() 함수를 활용하지만 Windows에서는 멀티 코어로 병렬 계산을 수행할 수 없다.
  • Windows 환경에서는 parApply() & parSapply() & parLapply() 함수를 활용한다.

parApply() & parSapply() & parLapply() 함수 사용법

• parApply() & parSapply() & parLapply() 함수 사용법에 앞서서 예시로 활용할 데이터는 다음과 같다.

##############################################
# library
##############################################
library(parallel)


##############################################
# data
##############################################
n <- 100
p <- 10000
Xmat <- matrix(rnorm(n*p), n, p)

• parApply() 함수의 사용법은 다음과 같다 :

• (Step 1) Cluster를 생성해야함. 먼저 사용하고 있는 컴퓨터 환경에서의 코어 개수를 먼저 확인할 것을 권장.

##############################################
# parApply
##############################################

### (Step 1) check the # of cores
numCores <- min( detectCores()-1, 12)

• (Step 2) Cluster를 등록해야함.

### (Step 2) set up a cluster 
myCluster <- makeCluster(numCores)

• (Step 3) 병렬 계산 수행.

### (Step 3) parApply
result <- parApply(myCluster, Xmat, 2, mean)

• (Step 4) 병렬 계산이 끝났으면 반드시 Cluster를 멈추어줄 것.

### (Step 4) stop the cluster 
stopCluster(myCluster)

• 이 때 한가지 또 주의할 점은 변수 스코프다.
  • parApply()을 사용할 때 변수 스코프 관련해서 접근할 수 없는 변수가 있는 경우가 잦음.
  • 특히 함수를 병렬 계산해서 사용할 때 변수 선언과 관련하여 이러한 문제가 잦다.
  • 이를 위해 clusterExport() 함수를 사용해주어야 한다.
  • 이는 cluster를 셋업했을 때 각 코어에서는 library나 object가 로드되어있지 않은 새로운 R세션에서 진행되기 때문에, export를 통해서 각 코어로 이들을 전달해주는 것이다.

##############################################
# function
##############################################
a <- 1
b <- 2
custoFUN0 <- function(x){ mean(x) + a*b }

##############################################
# parApply
##############################################

### (Step 1) check the # of cores
numCores <- min( detectCores()-1, 12)

### (Step 2) set up a cluster 
myCluster <- makeCluster(numCores)

### (Step 3-0) export to clusters
clusterExport(myCluster, varlist=c("a","b"))

### (Step 3) parApply
result <- parApply(myCluster, Xmat, 2, custoFUN0)

### (Step 4) stop the cluster 
stopCluster(myCluster)

• parSapply() 함수와 parLapply() 함수의 사용법도 parApply() 함수와 거의 동일하다. apply()대신 sapply()와 lapply() 함수 문법을 따라주면 된다.

##############################################
# parSapply
##############################################

### (Step 1) check the # of cores
numCores <- min( detectCores()-1, 12)

### (Step 2) set up a cluster 
myCluster <- makeCluster(numCores)

### (Step 3-0) export to clusters
clusterExport(myCluster, varlist=c("a","b","custoFUN0","Xmat"))

### (Step 3) parApply
result <- parSapply(myCluster, 1:ncol(Xmat), function(i){ custoFUN0(Xmat[,i]) })

### (Step 4) stop the cluster 
stopCluster(myCluster)


##############################################
# parLapply
##############################################

### (Step 1) check the # of cores
numCores <- min( detectCores()-1, 12)

### (Step 2) set up a cluster 
myCluster <- makeCluster(numCores)

### (Step 3-0) export to clusters
clusterExport(myCluster, varlist=c("a","b","custoFUN0","Xmat"))

### (Step 3) parApply
result <- parLapply(myCluster, 1:ncol(Xmat), function(i){ custoFUN0(Xmat[,i]) })

### (Step 4) stop the cluster 
stopCluster(myCluster)

apply() vs parApply()

• apply() 함수와 parApply() 함수의 성능을 비교해볼 것이다.

• (Case 1)

  • 다음에서 볼 수 있다시피 mean 함수 처럼 단순한 함수의 경우 apply()함수도 충분히 효율적으로 잘 작동하며 오히려 parApply()함수보다 더 좋은 성능을 나타낼 수 있음을 확인할 수 있다.

##############################################
# library
##############################################
library(parallel)
library(tictoc)


##############################################
# data
##############################################
n <- 100
p <- 10000
Xmat <- matrix(rnorm(n*p), n, p)


##############################################
# custom function
##############################################

### function 1
custoFUN1 <- function(x){ 
  
  # temporary setting
  tmp.n <- 1000
  tmp.p <- 5
  
  # temporary data
  tmp.y <- rnorm(tmp.n)
  tmp.x <- matrix(rnorm(tmp.n*tmp.p), tmp.n, tmp.p)
  
  # lm
  tmp.lm <- lm(tmp.y ~ tmp.x)
  tmp.result <- mean((fitted.values(tmp.lm) - tmp.y)^2)
  
  return(tmp.result)
}

### function 2
n0 <- 1000
p0 <- 5
custoFUN2 <- function(x){ 
  
  # temporary setting
  tmp.n <- n0
  tmp.p <- p0
  
  # temporary data
  tmp.y <- rnorm(tmp.n)
  tmp.x <- matrix(rnorm(tmp.n*tmp.p), tmp.n, tmp.p)
  
  # lm
  tmp.lm <- lm(tmp.y ~ tmp.x)
  tmp.result <- mean((fitted.values(tmp.lm) - tmp.y)^2)
  
  return(tmp.result)
}

##############################################
# apply vs parApply - case 1
##############################################

### apply 
tic("apply : ")
result0 <- apply(Xmat, 2, mean)
toc()

apply : : 0.05 sec elapsed

### parApply
tic("parApply : ")
# check the # of cores
tic("step 1 - ")
numCores <- min( detectCores()-1, 12)
toc()

step 1 - : 0 sec elapsed

# set up a cluster 
tic("step 2 - ")
myCluster <- makeCluster(numCores)
toc()

step 2 - : 0.22 sec elapsed

# parApply
tic("step 3 - ")
result <- parApply(myCluster, Xmat, 2, mean)
toc()

step 3 - : 0.04 sec elapsed

# stop the cluster 
tic("step 4 - ")
stopCluster(myCluster)
toc()

step 4 - : 0 sec elapsed

toc()

parApply : : 0.26 sec elapsed

• (Case 2)
  • 그러나 복잡한 함수를 활용할 경우 병렬계산의 이점이 나온다. 훨씬 더 빨리 계산됨을 확인할 수 있다.

##############################################
# apply vs parApply - case 2
##############################################

### apply 
tic("apply : ")
result0 <- apply(Xmat, 2, custoFUN1)
toc()

apply : : 6.05 sec elapsed

### parApply
tic("parApply : ")
# check the # of cores
tic("step 1 - ")
numCores <- min( detectCores()-1, 12)
toc()

step 1 - : 0 sec elapsed

# set up a cluster 
tic("step 2 - ")
myCluster <- makeCluster(numCores)
toc()

step 2 - : 0.2 sec elapsed

# parApply
tic("step 3 - ")
result <- parApply(myCluster, Xmat, 2, custoFUN1)
toc()

step 3 - : 0.75 sec elapsed

# stop the cluster 
tic("step 4 - ")
stopCluster(myCluster)
toc()

step 4 - : 0 sec elapsed

toc()

parApply : : 0.95 sec elapsed

sapply() vs parSapply()

##############################################
# sapply vs parSapply - case 2
##############################################

### sapply 
tic("sapply : ")
result0 <- sapply(1:ncol(Xmat), function(i){custoFUN2(Xmat[,i])})
toc()

sapply : : 5.45 sec elapsed

### parSapply
tic("parSapply : ")
# check the # of cores
tic("step 1 - ")
numCores <- min( detectCores()-1, 12)
toc()

step 1 - : 0 sec elapsed

# set up a cluster 
tic("step 2 - ")
myCluster <- makeCluster(numCores)
toc()

step 2 - : 0.21 sec elapsed

# export to clusters
tic("step 3-0 - ")
clusterExport(myCluster, varlist=c("n0","p0","custoFUN2"))
toc()

step 3-0 - : 0.01 sec elapsed

# parSapply
tic("step 3 - ")
result <- parSapply(myCluster, 1:ncol(Xmat), function(i){custoFUN2(Xmat[,i])})
toc()

step 3 - : 0.74 sec elapsed

# stop the cluster 
tic("step 4 - ")
stopCluster(myCluster)
toc()

step 4 - : 0 sec elapsed

toc()

parSapply : : 0.96 sec elapsed

lapply() vs parLapply()

##############################################
# sapply vs parLapply - case 2
##############################################

### lapply 
tic("lapply : ")
result0 <- lapply(1:ncol(Xmat), function(i){custoFUN2(Xmat[,i])})
toc()

lapply : : 5.65 sec elapsed

### parLapply
tic("parLapply : ")
# check the # of cores
tic("step 1 - ")
numCores <- min( detectCores()-1, 12)
toc()

step 1 - : 0 sec elapsed

# set up a cluster 
tic("step 2 - ")
myCluster <- makeCluster(numCores)
toc()

step 2 - : 0.22 sec elapsed

# export to clusters
tic("step 3-0 - ")
clusterExport(myCluster, varlist=c("n0","p0","custoFUN2"))
toc()

step 3-0 - : 0 sec elapsed

# parLapply
tic("step 3 - ")
result <- parLapply(myCluster, 1:ncol(Xmat), function(i){custoFUN2(Xmat[,i])})
toc()

step 3 - : 0.67 sec elapsed

# stop the cluster 
tic("step 4 - ")
stopCluster(myCluster)
toc()

step 4 - : 0 sec elapsed

toc()

parLapply : : 0.89 sec elapsed

2.2 foreach 패키지

• parApply() 관련 함수를 통하여 병렬 계산을 하는 것은 제한적인 부분이 있음

• 바로 apply 관련 함수의 문법을 따르기 때문에 복잡한 계산과정을 거치는 경우 적용하기에 불편한 부분이 있음.

• 대신에 for문을 작성하여 이를 병렬 계산하는 방법이 있음.

• 이를 위하여 foreach패키지를 활용하며, 이 때 병렬계산을 위하여 doParallel 패키지도 같이 활용. (parallel패키지를 foreach패키지와 함께 활용하기 위하여 변형된 패키지.)

• (예시 1)

  • %do%는 병렬화를 사용하지 않고 계산하지 않는 경우. 일반적인 for문과 같음
  • %dopar%는 병렬화를 사용하고 계산하는 경우.
    • 다만 이대로는 병렬화를 사용할 수 없고, doParallel패키지를 활용한 추가적인 작업이 필요.

######################################################
# library
######################################################
library(foreach)
library(doParallel)

######################################################
# foreach package
######################################################

### for loop
for(i in 1:3){ sqrt(i) }

### Without parallelization --> %do%
output <- foreach(i = 1:3) %do% {
  sqrt(i)
}
output

[[1]]
[1] 1

[[2]]
[1] 1.414214

[[3]]
[1] 1.732051

### With parallelization --> %dopar%
output <- foreach(i = 1:3) %dopar% {
  sqrt(i)
}

• foreach 문을 쓸 때 여러가지 option이 있는데 그중 하나가 .combine :

### foreach - other options : .combine
# .combine='c'
output <- foreach(i = 1:3, .combine='c') %do% {
  sqrt(i)
}
output

[1] 1.000000 1.414214 1.732051

# .combine='cbind' 
output <- foreach(x = 1:3, .combine = 'cbind') %do% {
  c(sqrt(x), log(x), x^2)
} 
output

     result.1  result.2 result.3
[1,]        1 1.4142136 1.732051
[2,]        0 0.6931472 1.098612
[3,]        1 4.0000000 9.000000

# .combine='rbind' 
output <- foreach(x = 1:3, .combine = 'rbind') %do% {
  c(sqrt(x), log(x), x^2)
} 
output

             [,1]      [,2] [,3]
result.1 1.000000 0.0000000    1
result.2 1.414214 0.6931472    4
result.3 1.732051 1.0986123    9

# .combine='+"
output <- foreach(i = 1:3, .combine='+') %do% {
  sqrt(i)
}
output

[1] 4.146264

# .combine='*"
output <- foreach(i = 1:3, .combine='*') %do% {
  sqrt(i)
}
output

[1] 2.44949

2.2 foreach 패키지 + doParallel 패키지

• 앞서서 설명했듯이, 병렬 계산을 하기 위해서는 %dopar%를 사용해야하는데, 이 때는 doParallel 패키지도 같이 활용해주어야 함.
  • 처음 병렬 계산을 셋업할 때 이 스텝에서 시간이 꽤 걸림
  • 따라서 시간이 얼마 걸리지 않는 for문의 경우, 이 때문에 오히려 병렬 계산이 시간이 더 걸릴 수 있음.
  • 그러나, 계산량이 많은 for문의 경우 병렬 계산이 압도적으로 효율적.
  • 단, R패키지 중에서 일부 함수는 병렬 계산이 수행되도록 코딩되어있는 경우도 있음. 해당 함수는 병렬 계산이 안될 것이기 때문에, R 패키지 활용 시 이 부분을 매뉴얼에서 확인해주어야할 필요성이 있음.
• (Step 1) Cluster를 생성해야함. 먼저 사용하고 있는 컴퓨터 환경에서의 코어 개수를 먼저 확인할 것을 권장.

######################################################
# foreach & doParallel packages
######################################################

### step 1
# check # of cores
detectCores()

[1] 32

# create the cluster
myCluster <- makeCluster(12)

• (Step 2) Cluster를 등록해야함.

### step 2
# register the cluster with the foreach package
registerDoParallel(myCluster)

• (Step 3) 병렬 계산 수행. 이 때, %do% 대신 %dopar% 사용

### step 3
# parallel computing using foreach with %dopar% 
output <- foreach(i = 1:3) %dopar% {
  sqrt(i)
}

• (Step 4) 병렬 계산이 끝났으면 반드시 Cluster를 멈추어줄 것.

### Step 4
# stop the cluster 
stopCluster(myCluster)

2.3 R코드 예시

######################################################
# R code exmaple 
######################################################
# foreach & doParallel packages
######################################################

tic("foreach & doParallel : ")
### step 1
# check # of cores
detectCores()

[1] 32

# create the cluster
myCluster <- makeCluster(12)

### step 2
# register the cluster with the foreach package
registerDoParallel(myCluster)

### step 3
# parallel computing using foreach with %dopar% 
output <- foreach(i = 1:200) %dopar% {
  
  n <- 100000
  p <- 20
  y <- sample(0:1, size=n, replace=T)
  x <- matrix(rnorm(n*p), nrow=n)
  glm.fit <- glm(y~x, family="binomial")
  glm.fit
}

### Step 4
# stop the cluster 
stopCluster(myCluster)
toc()

foreach & doParallel : : 25.41 sec elapsed

######################################################
# for loop
######################################################

tic("for loop : ")
output2 <- list()
for(i in 1:200){
  
  n <- 100000
  p <- 20
  y <- sample(0:1, size=n, replace=T)
  x <- matrix(rnorm(n*p), nrow=n)
  glm.fit <- glm(y~x, family="binomial")
  output2[[i]] <- glm.fit
}
toc()

for loop : : 82.66 sec elapsed

• 병렬 계산을 수행할 때, 다른 R패키지의 함수를 사용하고 싶을 때는, foreach의 옵션 중에서 .package로 해당 R패키지를 선언해주어야함
  • foreach 문 밖에서 해당 R 패키지를 불러왔더라도 foreach 문에서 다시 한번 선언해야함.

######################################################
# R code exmaple 
######################################################
# foreach & doParallel packages
######################################################

### step 1
# check # of cores
detectCores()

[1] 32

# create the cluster
myCluster <- makeCluster(3)

### step 2
# register the cluster with the foreach package
registerDoParallel(myCluster)

### step 3
# parallel computing using foreach with %dopar% 
output <- foreach(i = 1:10, .packages = c("glmnet")) %dopar% {
  
  n <- 10000
  p <- 10
  y <- rnorm(n)
  x <- matrix(rnorm(n*p), nrow=n)
  
  train <- sample(1:nrow(x), nrow(x)/2)
  test <- (-train)
  y.test <- y[test]
  
  ridge.mod <- glmnet(x[train,], y[train], alpha = 1, lambda = 1, thresh = 1e-12)
  ridge.pred <- predict(ridge.mod, newx=x[test, ])
  mean((ridge.pred-y.test)^2)
}

### Step 4
# stop the cluster 
stopCluster(myCluster)

• R 환경에서 병렬 계산을 할 때 변수 스코프 관련해서 foreach 문 안에서 접근할 수 없는 변수가 있는 경우가 잦음.
  • 특히 함수 안에서 병렬계산을 할 때 이런 경우가 많음
  • 이 때 .export에서 해당 변수를 선언해주면 이를 해결할 수 있음.

tmp.ParallelComputing___FUN <- function(x,y, niter, core.num=8){
  
  # create the cluster
  myCluster <- makeCluster(core.num)
  
  ### step 2
  # register the cluster with the foreach package
  registerDoParallel(myCluster)
  
  ### step 3
  # parallel computing using foreach with %dopar% 
  output <- foreach(i = 1:niter, .combine="c", .export=c("tmp.intercept1","tmp.intercept2")) %dopar% {
    
    tmp.result <- x + i*y + tmp.intercept1 * tmp.intercept2
  }
  
  ### Step 4
  # stop the cluster 
  stopCluster(myCluster)
  
  return(output)
}

tmp.intercept1 <- -1
tmp.intercept2 <- -3

tmp.result <- tmp.ParallelComputing___FUN(x=3, y=2, niter=30, core.num=3)
tmp.result

 [1]  8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50
[23] 52 54 56 58 60 62 64 66

Section 3. Rcpp package

• R의 느린 연산속도를 빠르게 하기 위한 다른 방법으로 C++을 이용해서 컴파일 하는 방법이 있음.
  • C++ 상에서 필요한 코드를 작성하고 다시 R로 불러들여서 더 빠른 연산을 수행하는 방식.
  • 이를 가능하게 하는 것이 Rcpp 패키지.
  • RcppArmadillo 및 RcppEigen 패키지를 같이 많이 활용함
    • Armadillo와 Eigen는 선형 대수 관련 C++ 라이브러리로서 이들을 R에서 손쉽게 사용하기 위해서 등장한 것이 RcppArmadillo 및 RcppEigen 패키지
    • https://arma.sourceforge.net/ 참조
    • https://eigen.tuxfamily.org/index.php?title=Main_Page 참ㅈ모
• Rcpp 패키지와 함께 병렬 계산을 수행하면 연산 속도를 비약적으로 상승시킬 수 있음.
• 튜토리얼 : https://statkim.github.io/rcpp-tutorial/ 참조.
  • https://lareale.tistory.com/303
  • http://adv-r.had.co.nz/Rcpp.html