【编程技巧（三）】并行化+分治：大矩阵运算的终极解决方案

在对一个大矩阵执行相关性计算或Jaccard Index的计算时，其实执行的是矩阵任意两行（这里假设要进行分析的对象是矩阵的每个行）之间的两两的计算，若这个矩阵的规模非常庞大，有n行时，计算的时间复杂度就是，这个时候可以采用并行化策略来加速这个进程（参考上文的 2. R中的并行化方法）：

StatOut <- parApply(cl, data, 1, fun, data)

这样就会实现将一个 n vs. n 的问题拆分成 n 个可以并行解决的 1 vs. n 的问题，当然通过设置线程数为，使得每次只并行执行m个 1 vs. n 的问题

然后再在函数fun内部再定义一个并行化计算来进一步并行化加速上面产生的 1 vs. n 的计算：

fun <- function(vec, data){
    ...
    parApply(cl, data, 1, fun2, vec)
}

在这个函数内部实现了将一个 1 vs. n 的问题拆分成 n 个可以并行解决的 1 vs. 1 的问题

这样就实现了两步并行化，这样在保证硬件条件满足的情况下，的确能显著加快分析速度

但是并行化技术会带来一个问题是：虽然时间开销减少了，但是空间开销显著增加了

比如，第一次并行化，将一个 n vs. n 的问题拆分成 次可以并行解决的 m个 1 vs. n 的问题，则需要在每一次并行化任务中拷贝一个 1 vs. n 的计算对象，即原始有n行的矩阵被拷贝了m次，则相应的缓存空间也增加了m倍，很显然内存的占用大大增加了

空间开销显著增加带来的后果就是，很容易导致运行内存不足程序运行中断的问题，那该怎么解决这个问题呢？

可以采用分治方法（Divide-Conquer)，将原始大矩阵，按照行拆分成多个小的子块，对每个子块执行计算，从而得到每个子块的运算结果，最后再讲每个子块的结果进行合并：

n_row <- nrow(data)
nblock <- 10 # 用于设定子块的数，用户可以自行指定
# 该列表变量用于保存每个子块的计算结果，若总共拆成了m个子块，
# 因为要进行任意两子块的计算，因此会有mxm个子块计算结果，因
# 此该列表要保存的数据为mxm维度的数据，每个元素为一个计算结
# 果，它们都是矩阵
StatOutList <- vector("list", nblock) 

# 1. 开始进行分块计算
print("[Divide-Conquer]Start carry out Divide-Conquer for Large Observed Matrix")
print("##################################################")
for(i in 1:nblock){
    for(j in 1:nblock){
        nrow_start <- (i-1)*nrow_block+1
        nrow_end <- i*nrow_block
        # 并行化计算
        if(!is.na(Args[2])){
            print(paste("[Divide-Conquer]Start carry out statistic Jaccard Index parallel for block: ",i,"-",j,sep=''))
            StatOutList[[i]] <- append(StatOutList[[i]], parApply(cl, data[nrow_start:nrow_end,], 1 , JaccardIndexSer, data))
            print(paste("[Divide-Conquer]Finish run parallel for block: ",i,"-",j,sep=''))
        # 串行计算
        }else{
            print(paste("[Divide-Conquer]Start carry out statistic Jaccard Index serially for block: ",i,"-",j,sep=''))
            StatOutList[[i]] <- append(StatOutList[[i]], apply(data, 1 , JaccardIndexSer, data))
            print(paste("[Divide-Conquer]Finish run serially for block: ",i,"-",j,sep=''))
        }
    }
}

# 2. 结束分治方法的分块计算
if(!is.na(Args[2])){
    print("##################################################")
    print("[Divide-Conquer]Finish parallel running for statistic Jaccard Index!")
    stopCluster(cl)
}else{
    print("##################################################")
    print("[Divide-Conquer]Finish serial running for statistic Jaccard Index!")
}

# 3. 开始进行子块结果的合并
print("[Divide-Conquer]Start bind sub-block statout")
StatOut <- vector("list", nblock)
# 先对列进行合并
for(i in 1:nblock){
    for(block in StatOutList[[i]]){
        StatOut[[i]] <- cbind(StatOut[[i]], block)
    }
}
# 再对行进行合并
StatOutMerge <- data.frame()
for(block in StatOut){
    StatOutMerge <- rbind(StatOutMerge, block)
}