R與并行計算

本文首先介紹了并行計算的基本概念，然后簡要闡述了R和并行計算的關(guān)系。之后作者從R語言用戶的使用角度討論了隱式和顯示兩種并行計算模式，并給出了相應(yīng)的案例。隱式并行計算模式不僅提供了簡單清晰的使用方法，而且很好的隱藏了并行計算的實現(xiàn)細節(jié)。因此用戶可以專注于問題本身。顯示并行計算模式則更加靈活多樣，用戶可以按照自己的實際問題來選擇數(shù)據(jù)分解，內(nèi)存管理和計算任務(wù)分配的方式。最后，作者探討了現(xiàn)階段R并行化的挑戰(zhàn)以及未來的發(fā)展。

R與并行計算

統(tǒng)計之都的小伙伴們對R，SAS，SPSS， MATLAB之類的統(tǒng)計軟件的使用定是輕車熟路了，但是對并行計算（又名高性能計算，分布式計算）等概念可能多少會感到有點陌生。應(yīng)太云兄之邀，在此給大家介紹一些關(guān)于并行計算的基本概念以及在R中的使用。

什么是并行計算？

并行計算，準確地說應(yīng)該包括高性能計算機和并行軟件兩個方面。在很長一段時間里，中國的高性能計算機處于世界領(lǐng)先水平。在最近一期的世界TOP500超級計算機排名中，中國的神威太湖之光列居榜首。 但是高性能計算機的應(yīng)用領(lǐng)域卻比較有限，主要集中在軍事，航天航空等國防軍工以及科研領(lǐng)域。對于多數(shù)個人，中小型企業(yè)來說高性能計算機還是陽春白雪。

不過，近年來隨著個人PC機，廉價機群，以及各種加速卡(NVIDIA GPU, Intel Xeon Phi, FPGA)的快速發(fā)展，現(xiàn)在個人電腦已經(jīng)完全可以和過去的高性能計算機相媲美了。相比于計算機硬件的迅速發(fā)展，并行軟件的發(fā)展多少有些滯后，試想你現(xiàn)在使用的哪些軟件是支持并行化運算的呢？ 

軟件的并行化需要更多的研發(fā)支持，以及對大量串行算法和現(xiàn)有軟件的并行化，這部分工作被稱之為代碼現(xiàn)代化（code modernization）。聽起來相當高大上的工作，然而在實際中大量的錯誤修正（BUGFIX），底層數(shù)據(jù)結(jié)構(gòu)重寫，軟件框架的更改，以及代碼并行化之后帶來的運行不確定性和跨平臺等問題極大地增加了軟件的開發(fā)維護成本和運行風險，這也使得這項工作在實際中并沒有想象中的那么吸引人。

R為什么需要并行計算？

那么言歸正傳，讓我們回到R本身。R作為當前最流行的統(tǒng)計軟件之一，具有非常多的優(yōu)點，比如豐富的統(tǒng)計模型與數(shù)據(jù)處理工具，以及強大的可視化能力。但隨著數(shù)據(jù)量的日漸增大，R的內(nèi)存使用方式和計算模式限制了R處理大規(guī)模數(shù)據(jù)的能力。從內(nèi)存角度來看，R采用的是內(nèi)存計算模式（In-Memory），被處理的數(shù)據(jù)需要預(yù)取到主存（RAM）中。其優(yōu)點是計算效率高、速度快，但缺點是這樣一來能處理的問題規(guī)模就非常有限（小于RAM的大?。Ａ硪环矫妫琑的核心（R core）是一個單線程的程序。因此，在現(xiàn)代的多核處理器上，R無法有效地利用所有的計算內(nèi)核。腦補一下，如果把R跑到具有260個計算核心的太湖之光CPU上，單線程的R程序最多只能利用到1/260的計算能力，而浪費了其他259/260的計算核心。

怎么破？并行計算！

并行計算技術(shù)正是為了在實際應(yīng)用中解決單機內(nèi)存容量和單核計算能力無法滿足計算需求的問題而提出的。因此，并行計算技術(shù)將非常有力地擴充R的使用范圍和場景。最新版本的R已經(jīng)將parallel包設(shè)為了默認安裝包?？梢奟核心開發(fā)組也對并行計算非常重視了。

R用戶：如何使用并行計算？

從用戶的使用方式來劃分，R中的并行計算模式大致可以分為隱式和顯示兩種。下面我將用具體實例給大家做一個簡單介紹。

隱式并行計算

隱式計算對用戶隱藏了大部分細節(jié)，用戶不需要知道具體數(shù)據(jù)分配方式 ，算法的實現(xiàn)或者底層的硬件資源分配。系統(tǒng)會根據(jù)當前的硬件資源來自動啟動計算核心。顯然，這種模式對于大多數(shù)用戶來說是最喜聞樂見的。我們可以在完全不改變原有計算模式以及代碼的情況下獲得更高的性能。常見的隱式并行方式包括：

1、 使用并行計算庫，如OpenBLAS，Intel MKL，NVIDIA cuBLAS

這類并行庫通常是由硬件制造商提供并基于對應(yīng)的硬件進行了深度優(yōu)化，其性能遠超R自帶的BLAS庫，所以建議在編譯R的時候選擇一個高性能庫或者在運行時通過LD_PRELOAD來指定加載庫。具體的編譯和加載方法可以參見這篇博客的附錄部分【1】。在下面左圖中的矩陣計算比較實驗中，并行庫在16核的CPU上輕松超過R原有庫百倍之多。在右圖中，我們可以看到GPU的數(shù)學(xué)庫對常見的一些分析算法也有相當顯著的提速。

2、使用R中的多線程函數(shù)

OpenMP是一種基于共享內(nèi)存的多線程庫，主要用于單節(jié)點上應(yīng)用程序加速。最新的R在編譯時就已經(jīng)打開了OpenMP選項，這意味著一些計算可以在多線程的模式下運行。比如R中的dist函數(shù)就 是一個多線程實現(xiàn)的函數(shù)，通過設(shè)置線程數(shù)目來使用當前機器上的多個計算核心，下面我們用一個簡單的例子來感受下并行計算的效率， GitHub上有完整代碼【2】， 此代碼需在Linux系統(tǒng)下運行。

#comparison of single thread and multiple threads run
for(i in 6:11) {
    ORDER <- 2^i
    m <- matrix(rnorm(ORDER*ORDER),ORDER,ORDER);
    .Internal(setMaxNumMathThreads(1)); .Internal(setNumMathThreads(1)); res <- system.time(d <- dist(m))
    print(res)
    .Internal(setMaxNumMathThreads(20)); .Internal(setNumMathThreads(20)); res <- system.time(d <- dist(m))
    print(res)
}

3、使用并行化包

在R高性能計算列表【3】中已經(jīng)列出了一些現(xiàn)有的并行化包和工具。用戶使用這些并行化包可以像使用其他所有R包一樣快捷方便，始終專注于所處理的問題本身，而不必考慮太多關(guān)于并行化實現(xiàn)以及性能提升的問題。我們以H2O.ai【4】為例。 H2O后端使用Java實現(xiàn)多線程以及多機計算，前端的R接口簡單清晰，用戶只需要在加載包之后初始化H2O的線程數(shù)即可，后續(xù)的計算， 如GBM，GLM, DeepLearning算法，將會自動被分配到多個線程以及多個CPU上。詳細函數(shù)可參見H2O文檔

【5】。

R is connected to the H2O cluster:
    H2O cluster uptime:         1 hours 53 minutes
    H2O cluster version:        3.8.3.3
    H2O cluster name:           H2O_started_from_R_patricz_ywj416
    H2O cluster total nodes:    1
    H2O cluster total memory:   1.55 GB
    H2O cluster total cores:    4
    H2O cluster allowed cores:  4
    H2O cluster healthy:        TRUE
    H2O Connection ip:          localhost
    H2O Connection port:        54321
    H2O Connection proxy:       NA
    R Version:                  R version 3.3.0 (2016-05-03)

顯示并行計算

顯式計算則要求用戶能夠自己處理算例中數(shù)據(jù)劃分，任務(wù)分配，計算以及最后的結(jié)果收集。因此，顯式計算模式對用戶的要求更高，用戶不僅需要理解自己的算法，還需要對并行計算和硬件有一定的理解。值得慶幸的是，現(xiàn)有R中的并行計算框架，如parallel (snow,multicores)，Rmpi和foreach等采用的是映射式并行模型（Mapping），使用方法簡單清晰，極大地簡化了編程復(fù)雜度。R用戶只需要將現(xiàn)有程序轉(zhuǎn)化為*apply或者for的循環(huán)形式之后，通過簡單的API替換來實現(xiàn)并行計算。對于更為復(fù)雜的計算模式，用戶可以通過重復(fù)映射收集（Map-Reduce）的過程來構(gòu)造。

下面我們用一元二次方程求解問題來介紹如何利用*apply和foreach做并行化計算，完整的代碼（ExplicitParallel.R）【6】可以在GitHuB上下載。首先，我們給出一個非向量化的一元二次方程求解函數(shù)，其中包括了對幾種特殊情況的處理，如二次項系數(shù)為零，二次項以及一次項系數(shù)都為零或者開根號數(shù)為負。我們隨機生成了3個大向量分別保存了方程的二次項，一次項和常數(shù)項系數(shù)。

# Not vectorized function
solve.quad.eq <- function(a, b, c) {
    # Not validate eqution: a and b are almost ZERO
    if(abs(a) < 1e-8 && abs(b) < 1e-8) return(c(NA, NA) )
    # Not quad equation
    if(abs(a) < 1e-8 && abs(b) > 1e-8) return(c(-c/b, NA))
    # No Solution
    if(b*b - 4*a*c < 0) return(c(NA,NA))
    # Return solutions
   x.delta <- sqrt(b*b - 4*a*c)
   x1 <- (-b + x.delta)/(2*a)
   x2 <- (-b - x.delta)/(2*a)
    return(c(x1, x2))
}

# Generate data
len <- 1e6
a <- runif(len, -10, 10)
a[sample(len, 100,replace=TRUE)] <- 0
b <- runif(len, -10, 10)
c <- runif(len, -10, 10)

apply實現(xiàn)方式： 首先我們來看串行代碼，下面的代碼利用lapply函數(shù)將方程求解函數(shù)solve.quad.eq映射到每一組輸入數(shù)據(jù)上，返回值保存到列表里。

# serial code
system.time(
    res1.s <- lapply(1:len, FUN = function(x) { solve.quad.eq(a[x], b[x], c[x])})
)

接下來，我們利用parallel包里的mcLapply （multicores）來并行化lapply中的計算。從API的接口來看，除了額外指定所需計算核心之外，mcLapply的使用方式和原有的lapply一致，這對用戶來說額外的開發(fā)成本很低。mcLapply函數(shù)利用Linux下fork機制來創(chuàng)建多個當前R進程的副本并將輸入索引分配到多個進程上，之后每個進程根據(jù)自己的索引進行計算，最后將其結(jié)果收集合并。在該例中我們指定了2個工作進程，一個進程計算1：(len/2), 另一個計算（len/2+1）：len的數(shù)據(jù)，最后當mcLapply返回時將兩部分結(jié)果合并到res1.p中。但是，由于multicores在底層使用了Linux進程創(chuàng)建機制，所以這個版本只能在Linux下執(zhí)行。

# parallel
library(parallel)
# multicores on Linux
system.time(
  res1.p <- mclapply(1:len, FUN = function(x) { solve.quad.eq(a[x], b[x], c[x])}, mc.cores = 2)
)

對于非Linux用戶來說，我們可以使用parallel包里的parLapply函數(shù)來實現(xiàn)并行化。parLapply函數(shù)支持Windows，Linux，Mac等不同的平臺，可移植性更好，但是使用稍微復(fù)雜一點。在使用parLapply函數(shù)之前，我們首先需要建立一個計算組（cluster）。計算組是一個軟件層次的概念，它指我們需要創(chuàng)建多少個R工作進程（parallel包會創(chuàng)建新的R工作進程，而非multicores里R父進程的副本）來進行計算，理論上計算組的大小并不受硬件環(huán)境的影響。比如說我們可以創(chuàng)建一個大小為1000的計算組，即有1000個R工作進程。 但在實際使用中，我們通常會使用和硬件計算資源相同數(shù)目的計算組，即每個R工作進程可以被單獨映射到一個計算內(nèi)核。如果計算群組的數(shù)目多于現(xiàn)有硬件資源，那么多個R工作進程將會共享現(xiàn)有的硬件資源。如下例我們先用detectCores確定當前電腦中的內(nèi)核數(shù)目。值得注意的是detectCores的默認返回數(shù)目是超線程數(shù)目而非真正物理內(nèi)核的數(shù)目。例如在我的筆記本電腦上有2個物理核心，而每個物理核心可以模擬兩個超線程，所以detectCores()的返回值是4。 對于很多計算密集型任務(wù)來說，超線程對性能沒有太大的幫助，所以使用logical=FALSE參數(shù)來獲得實際物理內(nèi)核的數(shù)目并創(chuàng)建一個相同數(shù)目的計算組。由于計算組中的進程是全新的R進程，所以在父進程中的數(shù)據(jù)和函數(shù)對子進程來說并不可見。因此，我們需要利用clusterExport把計算所需的數(shù)據(jù)和函數(shù)廣播給計算組里的所有進程。最后parLapply將計算平均分配給計算組里的所有R進程，然后收集合并結(jié)果。

#Cluster on Windows
cores <- detectCores(logical = FALSE)
cl <- makeCluster(cores)
clusterExport(cl, c('solve.quad.eq', 'a', 'b', 'c'))
system.time(
   res1.p <- parLapply(cl, 1:len, function(x) { solve.quad.eq(a[x], b[x], c[x]) })
)
stopCluster(cl)

for實現(xiàn)方式： for循環(huán)的計算和*apply形式基本類似。在如下的串行實現(xiàn)中，我們提前創(chuàng)建了矩陣用來保存計算結(jié)果，for循環(huán)內(nèi)部只需要逐一賦值即可。

# serial code
res2.s <- matrix(0, nrow=len, ncol = 2)
system.time(
    for(i in 1:len) {
        res2.s[i,] <- solve.quad.eq(a[i], b[i], c[i])
    }
)

對于for循環(huán)的并行化，我們可以使用foreach包里的%dopar% 操作將計算分配到多個計算核心。foreach包提供了一個軟件層的數(shù)據(jù)映射方法，但不包括計算組的建立。因此，我們需要doParallel或者doMC包來創(chuàng)建計算組。計算組的創(chuàng)建和之前基本一樣，當計算組建立之后，我們需要使用registerDoParallel來設(shè)定foreach后端的計算方式。這里我們從數(shù)據(jù)分配方式入手，我們希望給每個R工作進程分配一段連續(xù)的計算任務(wù)，即將1：len的數(shù)據(jù)均勻分配給每個R工作進程。假設(shè)我們有兩個工作進程，那么進程1處理1到len/2的數(shù)據(jù)，進程2處理len/2+1到len的數(shù)據(jù)。所以在下面的程序中，我們將向量均勻分配到了計算組，每個進程計算chunk.size大小的聯(lián)系任務(wù)。并且在進程內(nèi)創(chuàng)建了矩陣來保存結(jié)果，最終foreach函數(shù)根據(jù).combine指的的rbind函數(shù)將結(jié)果合并。

# foreach
library(foreach)
library(doParallel)
# Real physical cores in the computer
cores <- detectCores(logical=F)
cl <- makeCluster(cores)
registerDoParallel(cl, cores=cores)
# split data by ourselves
chunk.size <- len/cores
system.time(
  res2.p <- foreach(i=1:cores, .combine='rbind') %dopar%
  {  # local data for results
     res <- matrix(0, nrow=chunk.size, ncol=2)
     for(x in ((i-1)*chunk.size+1):(i*chunk.size)) {
        res[x - (i-1)*chunk.size,] <- solve.quad.eq(a[x], b[x], c[x])
     }
     # return local results
     res
  }
)
stopImplicitCluster()
stopCluster(cl)

最后，我們在Linux平臺下使用4個線程進行測試，以上幾個版本的并行實現(xiàn)均可達到3倍以上的加速比。

R并行化的挑戰(zhàn)與展望

挑戰(zhàn)：

在實際中，并行計算的問題并沒有這么簡單。要并行化R以及整個生態(tài)環(huán)境的挑戰(zhàn)仍然巨大。

1、R是一個分散的，非商業(yè)化的軟件

R并不是由一個緊湊的組織或者公司開發(fā)的，其大部分包是由用戶自己開發(fā)的。這就意味著很難在軟件架構(gòu)和設(shè)計上來統(tǒng)一調(diào)整和部署。一些商業(yè)軟件，比如Matlab，它的管理維護開發(fā)就很統(tǒng)一，架構(gòu)的調(diào)整和局部重構(gòu)相對要容易一些。通過幾個版本的迭代，軟件整體并行化的程度就要高很多。

2、R的底層設(shè)計仍是單線程，上層應(yīng)用包依賴性很強

R最初是以單線程模式來設(shè)計的，意味著許多基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)并不是線程安全的。所以，在上層并行算法實現(xiàn)時，很多數(shù)據(jù)結(jié)構(gòu)需要重寫或者調(diào)整，這也將破壞R本來的一些設(shè)計模式。另一方面，在R中，包的依賴性很強，我們假設(shè)使用B包，B包調(diào)用了A包。如果B包首先實現(xiàn)了多線程，但是在一段時間之后A包也做了并行化。那么這時候就很有可能出現(xiàn)混合并行的情況，程序就非常有可能出現(xiàn)各種奇怪的錯誤（BUG），性能也會大幅度降低。

展望：

未來 R并行化的主要模式是什么樣的？以下純屬虛構(gòu)，如果有雷同完全是巧合。

1、R將會更多的依賴于商業(yè)化以及研究機構(gòu)提供的高性能組件

比如H2O，MXNet，Intel DAAL，這些包都很大程度的利用了并行性帶來的效率提升，而且有相關(guān)人員長期更新和調(diào)優(yōu)。從本質(zhì)上來說，軟件的發(fā)展是離不開人力和資金的投入的。

2、云計算平臺

隨著云計算的興起，數(shù)據(jù)分析即服務(wù)（DAAS：Data Analyst as a Services）以及機器學(xué)習即服務(wù)（MLAS： machine learning as a services）的浪潮將會到來。 各大服務(wù)商從底層的硬件部署，數(shù)據(jù)庫優(yōu)化到上次的算法優(yōu)化都提供了相應(yīng)的并行化措施，比如微軟近期推出了一系列R在云上的產(chǎn)品，更多信息請參見這篇文章。因此，未來更多的并行化工作將會對用戶透明，R用戶看到的還是原來的R，然而真正的計算已經(jīng)分布到云端了。