用R语言实现向量化与并行计算

应用场景决定知识的储备与工具的选择，反过来，无论你选择了什么样的工具，你一定会努力地把它改造成符合自己应用场景所需的那个样子。从这个道理来说，我选择了R作为数据挖掘人员手中攻城陷池的那把云梯，并努力地把它改造成自己希望的那个样子。

我最初接触到专门用于科学计算的工具，是大名鼎鼎的matlab，正如它帮助了无数中国学生顺利毕业的赫赫功劳一样，它是我对于向量化计算的启蒙老师。用过matlab的人都会对其循环结构的效率无法忍受，不知道是否有意而为之的这样的设计缺陷，迫使人们要想真正地用好它，就得接受它提供的向量化计算的思想，在掌握了这个专门为高效计算而设计的计算思想之后，你会发现自己获得的不单是计算效率上的极大提升，还有算法设计思想上高屋建瓴式的跨越。

毕竟matlab是更适合于研究领域的商业软件，而且是闭源的，毕业后不久，我就选择了R作为matlab的替代品，看中的正是R在向量化计算支持之余的灵活性及丰富的第三方库，似乎天生就是数据挖掘人员最趁手的那把刀。因为我需要的就是这样的一个计算环境，它不单是一门编程语言，也不必是一个已经很完备的工具。它就是这样的一个环境，拥有很多的各领域相关的工具包，我可以不必操心太多过于底层的或与工作主题没有直接关系的问题，同时当不满意时我又具有对它最自由的掌控。实际上，我寻求的就是matlab与C之间的一个平衡点。R是一个面向科学工程计算特别是统计计算的工具，与matlab一样，其循环结构的效率也无法让人满意，所以，我们必须利用向量化的编程范式，必要时采用并行/分布计算（因为，向量化本质上就是一种并行计算，也是我们通常理解那种并行计算的天然先驱）。而这一切，在R面前，都不是什么问题。

所谓向量化：

wikipedia中的定义是：Vectorization is the more limited process of converting a computer program from a scalar implementation, which processes a single pair of operands at a time, to a vector implementation which processes one operation on multiple pairs of operands at once。向量化计算是一种特殊的并行计算的方式，相比于一般程序在同一时间只执行一个操作的方式，它可以在同一时间执行多次操作，通常是对不同的数据执行同样的一个或一批指令，或者说把指令应用于一个数组/向量。

像R和matlab这样的现代科学软件包，都会以向量化作为其计算的基本特点（即使python的numpy包也是如此），所以在R的基本运算中，随处可见向量化计算的影子，以下仅举几例，以让读者了解向量化是多么地简单和直观：

向量取值：V[1:10]  （把get操作应用于向量V的不同元素）

向量赋值：V[1:10] <- seq(1,10)  （把set操作应用于一个序列与向量V的对应元素）

apply系列：lapply(V, mean)  （跟python的map函数类似，是向量化最直接的表达形式）

矩阵运算：A + B；A %*% B  （矩阵的基本运算也是向量化的典型形式）

这些操作很常见，以致于我们都没有意识到这就是一种有助于提高计算性能的向量化计算，更忘了由此而把这样的思想扩展到我们算法设计的更多方面。熟练使用它，你获得不单是计算上的好处，还有对问题理解上的整体性。

向量化的处理方式是现代计算机的一个特点，无论硬件还是软件上，都提供了支持。

硬件上的支持：Intel’s MMX, SSE, and ARM’s NEON instruction sets support such vectorized loops.（摘自wikipedia：http://en.wikipedia.org/wiki /Vectorization_(computer_science)）

软件上的支持：著名的线性代数运算包blas对矩阵运算的自动并行化。例如，在R里执行如下的简单命令，在不需要作任何额外工作的情况下，就可以自动地实现并行化（前提是你的机器安装了blas）。

 M = matrix(5000*5000,5000,5000) ## 生成一个5000*5000的矩阵
S = M %*% M ## 作矩阵乘法，在多核并安装了blas的机器里，会自动并行

应用向量化的思维方式去解决问题：

当我们习惯了用C语言的单步循环的思想来思考问题的时候，要把思维切换到向量化的方式是件比较困难的事件，以下显示一个例子，看用向量化的思想是怎么去解决问题，这样的描述又是多么美。

任务：对于稀疏矩阵M，让其第i(i＝1…m)行的各非零元素减去某个值w[i]。

正常的循环式思维的解决方案比较直观，作两层循环，让第i行非零元素减去w[i]即可，但这样的操作在脚本语言特别是像R这样的科学计算语言里执行效率不高，而当你面对的是一个稀疏矩阵时，问题又会变得复杂。

向量化的解决方案如下：

 ## 利用向量化的计算范式设计算法
## 实例：稀疏矩阵不同行的各元素减去某个值  
library(Matrix)
generate <- function(nr, nc, sr){
L = nr*nc
M = Matrix(data=0, nrow=nr, ncol=nc)
idx = sample(1:L,as.integer(sr*L))
val = rnorm(sr*L, mean=1)
M[idx] = val
return(M)
}

M = generate(10000,5000,0.1)  ## 生成稀疏矩阵
V = rnorm(10000)  ## 矩阵各行减去的值
N = M
N@x = 1  ## 复制并设置矩阵
A = Diagonal(nrow(N), V)%*% N  ## 把V的值放到新矩阵的各行非零单元中
M@x = M@x- A@x  # M= M-A

除去数据准备阶段，真正用于解决问题的只有第16到19行的4行命令，十分短小精悍，但需要一定的线性代数知识去理解这个过程。

从向量化到并行计算：

除了计算机硬件与科学计算包的支持外，向量化计算还是并行计算的天然先驱，如果你向量化后的算法效率仍然不佳，可以考虑把它并行化。

R那浩如烟海的第三方包里提供了大量支持并行计算的包，这里选择的是snowfall这个包，它可以简单地构建一个计算集群，把计算并行分布到集群的不同节点进行计算。以下通过一个例子比较循环，向量化以及并行三种方式在效率上的差异。

实战案例：for循环，apply，snowfall（并行）的比较。

任务：对一个矩阵每列排序并取回前10个。

 library(snowfall) 
mysort <- function(x){
  # replicate(2, sort(x))
   return(sort(x)[1:10])
}  

do_for <- function(M){
  ans = matrix(0,10, ncol(M))
  for(iin1:ncol(M)){
    ans[,i]= sort(M)[1:10]
  }
}  

do_apply <- function(M){
  return(apply(M,2, mysort))
}  

do_snowfall <- function(M, ncl){
  sfInit(parallel=TRUE, cpus=ncl) ##初始化集群环境
  ans <- sfApply(M,2, mysort) ##把任务分发到各个slave
  sfStop() ##关闭集群
  return(ans)
}  

#M = matrix(rnorm(10000000),100,100000)
#system.time(ans<- do_for(M))
#system.time(ans<- do_apply(M))
#system.time(ans<- do_snowfall(M,4)) ##用4个核并行

在我的服务器里运行的结果是这样的，do_for循环基本不可用，do_apply可以在一分钟内运算完毕，而do_snowfall的时间为 do_apply的一半。当取消第4行的注释，即增加各个子任务的计算负荷后，do_apply与do_snowfall的计算时间都增加，而 do_snowfall的时间变为do_apply的三分之一。

结论1：向量化的计算方式比起传统的循环计算有极大的性能优势。

结论2：由于并行的过程为：master把任务分解，分发到多个slave进行运算，slave返回结果到master。所以多核计算并不一定比最优的单核计算快，要看性能的瓶颈在slave还是在master上。

结论3：适合并行/分布计算的情景主要有两种：一是各slave的计算耗费为瓶颈，并行到多核能减少计算时间，越是slave耗时型的计算并行收益越大；二是一台机器的内存不足以进行整体的计算，分布到多台机器计算能把内存占用分开，这种情况下即使分布计算比单机慢也是可以接受的。

与map-reduce扯扯关系：

Map-reduce是google三驾马车之一，试图把所有计算都转换为map和reduce两种基本的运算方式，而达到良好的可并行性，从而实现计算上的可扩展性。虽然并不是每个公司都能实现像map-reduce那样的巨型计算框架，但是如果你能熟练地使用向量化的思想设计你的算法，那其实就是一个map-reduce的超集。向量化计算这个编程范式比map-reduce要复杂，但本质上都是使用了两种基本的运算，如果你把apply想像为 map，把矩阵或矢量的乘运算想像为reduce，它们就是一体的。你总可以把向量化的计算通过apply和矩阵运算来实现，如果有一天你拥有了一个类 map-reduce的计算框架，根据它们的对应关系，你的算法迁移就会非常地方便，甚至，你可以在你已经实现的向量化算法集当中抽取出一些共性来搭建自己的map-reduce系统。

别忘了一切优化的终极准则：过早的优化是魔鬼。

并行计算或分布式计算是为数据量与计算量日益膨胀的情况下所必须考虑的，但它的逻辑结构与维护成本也要高得多，如果你的应用场景还没达到需要进行并行或者分布计算的程度，没必要淌这趟混水。但是，向量化计算的编程范式却是很重要的，一个用向量化计算方式实现的算法，不但在当时获得了效益，而且其可扩展性也必定是强悍的，因为正如前文所述，向量化计算就是并行计算的天然先驱。除此之外，经常用向量化的方式来思考你的问题，你一定能感受到一种整体之美。

后记：

以上内容是从我参加第三届中国R语言会议的讲稿中整理出来的，以讲向量化的思想为主，R语言为辅，当时的PPT可以从这里下载（PPT 上内容比较少，可以下载下来看备注部分）。我一直认为，通过这样的讲演，除了布道，更重要的是能够对自己一段时间以来的思考与工作做一个总结，其间无论是从自身整理还是从他人的提问，自己从中都能有所收获。根据当场及事后一些听众的疑问，我后续整理了下面一些内容。

如果问题很难向量化怎么办：并不是说R中不能使用for式的循环，事实上只有一层的for循环，并且对性能没有太大的影响的话，是不需要作过多的优化的。如果是两层for循环，可以考虑用apply去掉一层，如果有三层的循环，那很有可能就是你的算法有问题了。如果你的问题确实很难转换成向量化的形式，又或者说你已经懒于这样去做，那么不妨把最难于转换的一部分用C来实现，然后给R做一个接口，剩下的再考虑向量化会简单很多。我提倡的是结合R自有函数与向量化思想进行编程，因为R的自有函数大多都是用底层语言实现好的，效率有保证。

有人疑问是不是所有算法都可以并行：非基本的算法或多或少都是可以并行的，像kmeans，虽然从整体来说是一个迭代依赖的非并行式算法，但在每一步的迭代中却是可以实施并行的。应用map-reduce框架基于这么一个假设，你的问题，总可以找到一种算法来实现，这种算法可以或大部分可以用map 加reduce基本操作来实现。如果向量化能跟map-reduce对应起来，那么也可以认为，你的问题，总能找到一种可以用向量化思想表示的算法来完全或部分解决，如果有些部分很难向量化，参考上一个问题。

我选择的算法本身就很复杂：有这么一个说法，如果你的数据不足够多，信息不够完整，你会想出很多奇思怪招来从中获得结论，实际上很多科学上精巧而复杂的算法是基于这种情况而设计出来的。而在实际中如果你的数据急剧膨胀，信息与噪声都足够地多，你的问题焦点就变为如何用一个或一些简单有效算法或算法的组合来提取有价值的信息。如果你在小数据集的时候习惯了采用复杂的算法，当你的应用场景转变了，数据规模变了，却仍旧沿用原有的思维方式，是注定要尴尬的（仅针对工程应用环境，而非科学应用环境）。

原文链接：http://www.wentrue.net/blog/?p=945