【R Special】Copy-on-Modify机制及R性能优化第一步-CSDN博客

周日应邀在第六届R会议做了一个R工程实践的topic，聊了聊自己在实际工作中使用R碰到的坑。由于我的讲稿一向很简单，怕别人拿到后看起来吃力。而Copy-on-Modify又是万坑之源，一种原罪级别的存在，感觉有必要拿出来单独说说，以帮助R的新手了解自己的代码为什么慢，以及解决这个问题的第一步方案是什么。之所以只说第一步，是因为这一步是普适的，也是能给人带来极大思维乐趣的。后续的优化步骤因案例而变，有人会选择用C/C++接口优化原有模块，有人选择并行化，又会碰到很多新的目前也未能解决得很好的问题，我们也正在这一步的摸索途中。

关于R的Copy-on-Modify机制，Hadley在【这里】有比较详细的叙述，用一个案例来说明，那就是：

> address <- function(x) .Internal(inspect(x))   ## 包装一个获取变量内存地址的函数以方便使用
> x <- 1:10
> address(x)  ## 向量x的内存地址
@154baf8 13 INTSXP g0c4 [NAM(2)] (len=10, tl=0) 1,2,3,4,5,...
> x[1] <- 2   ## 改变了x中的某个值
> address(x)  ## x的内存地址发生了拷贝
@d10ab8 14 REALSXP g0c6 [NAM(2)] (len=10, tl=0) 2,2,3,4,5,...

注：以上例子由于没有考虑到数据类型的原因，已经不适用，详情请看comment里的讨论。而且随着R版本的迭代，copy-on-modify的情况可能会越来越少，让我们共同期待这样的改进。

非常显而易见的机制：当需要对变量x自身或x中某个单元赋值，R会把x拷贝到新的内存地址，然后赋值。这种机制与很多流行工程语言中的对象引用机制不一样，比如如下的python代码：

>>> x = range(10)
>>> y = x
>>> y
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> x[1] = 2
>>> y
[0, 2, 2, 3, 4, 5, 6, 7, 8, 9]

x和y指向同一块内存区域，x内的单元发生改变时，是直接对这块内存的操作，所以x和y同时改变。这种情况发生在R的身上，就是y仍然指向原来的内存块，x被拷贝后重赋值，两者从内存到取值都已经不同（注：x的拷贝只发生在x被修改时，即on modify）。

比较明显的是，copy-on-modify机制带来了安全，你不用担心多个变量引用同一块内存区域的情况，以及出现类似于指针的evil。但它同时带来了极大的性能损耗——每一次变量修改时，都会产生一次内存拷贝。如果这个变量是一个大规模的矩阵，这种损耗带来的影响将是巨大的。所以R性能优化的第一步，也是必须要做的一步，就是尽量避免循环，把循环转换成向量化计算，因为向量化是一个整体的计算，计算前后只发生一次内存拷贝。R像其他科学计算语言一样原生提供了对向量化的支持，最基本的取值赋值都支持向量化的操作方式，再比如apply系列的函数，各种矩阵运算。关于R的向量化，在我以前的文章有过论述。

下面举一个例子来说明把循环逻辑向量化的必要性，例子取材于我面试时常用的一个题目：稀疏矩阵M的每一列非零元素减去该列的均值。

最直接的想法是这样的：

cs = colMeans(M)
for (i in nrow(M))
    for (j in ncol(M))
        M[i, j]=M[i, j]-cs[j]

我们算一下，根据copy-on-modify的原理，M矩阵在这个计算过程中被拷贝了nrow(M) * ncol(M)次，这还没算频繁拷贝导致的内存清理（gc）带来的额外损耗。你说这样的程序能不慢么。好，接下来用向量化的思维去改进一下。

cs = colMeans(M)
for (i in ncol(M))
    M[i, M[i,]!=0]=M[i, M[i,]!=0]-cs

这样拷贝的次数降为ncol(M)次，还是很多。而且考虑到这是一个稀疏矩阵，有特殊的存储结构，把它内部的作为个体去对待，单个取值赋值本身也是个性能大户，而矩阵的运算则是每一个稀疏矩阵库都必然会提供高效支持的，于是考虑把这个问题矩阵化，用矩阵运算的思维去解决。

library(Matrix)
...
N=M
## copy-on-modify，把M的数据拷贝一份，由N指向，所有非零元赋值为1
N@x = 1
## 后半部分是得到均值向量，并由此生成对角矩阵，最后右乘N的结果有兴趣的可以自行推演
S = N %*% Diagonal(colSums(M)/colSums(N))
M = M-S

这里使用了R的Matrix包，它对稀疏数据的存储结构可参考【这里】的csc格式。这个算法对M的拷贝次数为两次。理解这个过程需要初步的线性代数知识，但理解起来并不难。这个解决方案中没有使用任何判断和循环的控制结构，也没有细致到去处理每一个单元的赋值，而是从整体的角度来思考和解决问题。这也是数学思维的一部分，让你把单纯的敲击键盘行为，变为有趣的思维过程，这是最初R及向量化计算吸引我的原因。

会议上我另一个议题是豆瓣的新角色：Data Scientist。有兴趣的可以看看flycondor的【日记】，那里有比较详细的介绍。

不过他那个日记太累赘了，能细细看完的那都是绝顶的真爱。他后来给我总结了三个工作内容，对应于讲稿中的三个案例：管理数据、回答问题、探索数据。我扩充一下，第一个任务是常规性工作，产出报表和描述性统计数据，问题定义清晰，实现路径也清晰；第二个任务跟具体产品相关，通过数学建模解决产品线关心的问题，问题定义清晰，实现路径不清晰；第三个任务属于由你的发现向产品线的推动，需要较高的数据素质，问题定义不清晰，实现路径也不清晰。三个工作对人的要求依次递进。