使用 C 和并行化在 R 中快速关联
Posted
技术标签:
【中文标题】使用 C 和并行化在 R 中快速关联【英文标题】：Fast correlation in R using C and parallelization 【发布时间】：2013-09-28 16:17:57 【问题描述】：我今天的项目是使用我拥有的基本技能在 R 中编写一个快速关联例程。我必须找到近 400 个变量之间的相关性，每个变量都有近一百万个观察值（即大小为 p=1MM 行和 n=400 列的矩阵）。
对于 1MM 行和每个变量 200 个观察值，R 的原生相关函数需要将近 2 分钟。我没有为每列运行 400 次观察，但我猜这需要将近 8 分钟。我还有不到 30 秒的时间来完成它。
因此，我想做一些事情。
1 - 用 C 语言编写一个简单的相关函数并将其并行应用到块中（见下文）。
2 - 块 - 将相关矩阵分成三个块（大小为 K*K 的左上正方形、大小为 (pK)(pK) 的右下正方形和大小为 K 的右上矩形矩阵 em>(pK))。这涵盖了相关矩阵corr 中的所有单元格，因为我只需要上三角。
3 - 通过 .C 调用使用降雪并行运行 C 函数。
n = 100
p = 10
X = matrix(rnorm(n*p), nrow=n, ncol=p)
corr = matrix(0, nrow=p, ncol=p)
# calculation of column-wise mean and sd to pass to corr function
mu = colMeans(X)
sd = sapply(1:dim(X)[2], function(x) sd(X[,x]))
# setting up submatrix row and column ranges
K = as.integer(p/2)
RowRange = list()
ColRange = list()
RowRange[[1]] = c(0, K)
ColRange[[1]] = c(0, K)
RowRange[[2]] = c(0, K)
ColRange[[2]] = c(K, p+1)
RowRange[[3]] = c(K, p+1)
ColRange[[3]] = c(K, p+1)
# METHOD 1. NOT PARALLEL
########################
# function to calculate correlation on submatrices
BigCorr <- function(x)
Rows = RowRange[[x]]
Cols = ColRange[[x]]
return(.C("rCorrelationWrapper2", as.matrix(X), as.integer(dim(X)),
as.double(mu), as.double(sd),
as.integer(Rows), as.integer(Cols),
as.matrix(corr)))
res = list()
for(i in 1:3)
res[[i]] = BigCorr(i)
# METHOD 2
########################
BigCorr <- function(x)
Rows = RowRange[[x]]
Cols = ColRange[[x]]
dyn.load("./rCorrelation.so")
return(.C("rCorrelationWrapper2", as.matrix(X), as.integer(dim(X)),
as.double(mu), as.double(sd),
as.integer(Rows), as.integer(Cols),
as.matrix(corr)))
# parallelization setup
NUM_CPU = 4
library('snowfall')
sfSetMaxCPUs() # maximum cpu processing
sfInit(parallel=TRUE,cpus=NUM_CPU) # init parallel procs
sfExport("X", "RowRange", "ColRange", "sd", "mu", "corr")
res = sfLapply(1:3, BigCorr)
sfStop()
这是我的问题：
对于方法 1，它有效，但不是我想要的方式。我相信，当我传递 corr 矩阵时，我传递了一个地址，而 C 将在源代码处进行更改。
# Output of METHOD 1
> res[[1]][[7]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0.1040506 -0.01003125 0.23716384 -0.088246793 0 0 0 0 0
[2,] 0 1.0000000 -0.09795989 0.11274508 0.025754150 0 0 0 0 0
[3,] 0 0.0000000 1.00000000 0.09221441 0.052923520 0 0 0 0 0
[4,] 0 0.0000000 0.00000000 1.00000000 -0.000449975 0 0 0 0 0
[5,] 0 0.0000000 0.00000000 0.00000000 1.000000000 0 0 0 0 0
[6,] 0 0.0000000 0.00000000 0.00000000 0.000000000 0 0 0 0 0
[7,] 0 0.0000000 0.00000000 0.00000000 0.000000000 0 0 0 0 0
[8,] 0 0.0000000 0.00000000 0.00000000 0.000000000 0 0 0 0 0
[9,] 0 0.0000000 0.00000000 0.00000000 0.000000000 0 0 0 0 0
[10,] 0 0.0000000 0.00000000 0.00000000 0.000000000 0 0 0 0 0
> res[[2]][[7]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 -0.02261175 -0.23398448 -0.02382690 -0.1447913 -0.09668318
[2,] 0 0 0 0 0 -0.03439707 0.04580888 0.13229376 0.1354754 -0.03376527
[3,] 0 0 0 0 0 0.10360907 -0.05490361 -0.01237932 -0.1657041 0.08123683
[4,] 0 0 0 0 0 0.18259522 -0.23849323 -0.15928474 0.1648969 -0.05005328
[5,] 0 0 0 0 0 -0.01012952 -0.03482429 0.14680301 -0.1112500 0.02801333
[6,] 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000
[7,] 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000
[8,] 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000
[9,] 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000
[10,] 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.0000000 0.00000000
> res[[3]][[7]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000
[3,] 0 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000
[4,] 0 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000
[5,] 0 0 0 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000
[6,] 0 0 0 0 0 1 0.03234195 -0.03488812 -0.18570151 0.14064640
[7,] 0 0 0 0 0 0 1.00000000 0.03449697 -0.06765511 -0.15057244
[8,] 0 0 0 0 0 0 0.00000000 1.00000000 -0.03426464 0.10030619
[9,] 0 0 0 0 0 0 0.00000000 0.00000000 1.00000000 -0.08720512
[10,] 0 0 0 0 0 0 0.00000000 0.00000000 0.00000000 1.00000000
但原来的corr矩阵保持不变：
> corr
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 0
问题 #1：有什么方法可以确保 C 函数在源代码中更改 corr 的值？我仍然可以合并这三个以创建上三角相关矩阵，但我想知道是否可以从源头进行更改。注意：这并不能帮助我实现快速关联，因为我只是在运行一个循环。
问题 #2：对于方法 2，我如何将共享对象加载到每个核心，以便在初始化步骤中每个核心上的并行作业（而不是我是如何做到的）？ p>
问题 #3：这个错误是什么意思？我需要一些指示，我很想自己调试它。
问题 #4：有没有一种快速的方法可以在 30 秒内计算 1MM x 400 矩阵的相关性？
当我运行 METHOD 2 时，我收到以下错误：
R(6107) malloc: *** error for object 0x100664df8: incorrect checksum for freed object - object was probably modified after being freed.
*** set a breakpoint in malloc_error_break to debug
Error in unserialize(node$con) : error reading from connection
下面附上我用于关联的普通 C 代码：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stddef.h>
#include <R.h> // to show errors in R
double calcMean (double *x, int n);
double calcStdev (double *x, double mu, int n);
double calcCov(double *x, double *y, int n, double xmu, double ymu);
void rCorrelationWrapper2 ( double *X, int *dim, double *mu, double *sd, int *RowRange, int *ColRange, double *corr)
int i, j, n = dim[0], p = dim[1];
int RowStart = RowRange[0], RowEnd = RowRange[1], ColStart = ColRange[0], ColEnd = ColRange[1];
double xyCov;
Rprintf("\n p: %d, %d <= row < %d, %d <= col < %d", p, RowStart, RowEnd, ColStart, ColEnd);
if(RowStart==ColStart && RowEnd==ColEnd)
for(i=RowStart; i<RowEnd; i++)
for(j=i; j<ColEnd; j++)
Rprintf("\n i: %d, j: %d, p: %d", i, j, p);
xyCov = calcCov(X + i*n, X + j*n, n, mu[i], mu[j]);
*(corr + j*p + i) = xyCov/(sd[i]*sd[j]);
else
for(i=RowStart; i<RowEnd; i++)
for (j=ColStart; j<ColEnd; j++)
xyCov = calcCov(X + i*n, X + j*n, n, mu[i], mu[j]);
*(corr + j*p + i) = xyCov/(sd[i]*sd[j]);
// function to calculate mean
double calcMean (double *x, int n)
double s = 0;
int i;
for(i=0; i<n; i++)
s = s + *(x+i);
return(s/n);
// function to calculate standard devation
double calcStdev (double *x, double mu, int n)
double t, sd = 0;
int i;
for (i=0; i<n; i++)
t = *(x + i) - mu;
sd = sd + t*t;
return(sqrt(sd/(n-1)));
// function to calculate covariance
double calcCov(double *x, double *y, int n, double xmu, double ymu)
double s = 0;
int i;
for(i=0; i<n; i++)
s = s + (*(x+i)-xmu)*(*(y+i)-ymu);
return(s/(n-1));
【问题讨论】：
@MartinMorgan - R 的原生 cor 函数（基于我拥有的构建）需要更多时间，正如我上面提到的。我在下面使用 Andrey 的建议，1MM x 400 vars 大约需要 2 分钟。会更新。 【参考方案1】：使用快速 BLAS（通过 Revolution R 或 Goto BLAS），您可以在 R 中快速计算所有这些相关性，而无需编写任何 C 代码。 在我的第一代 Intel i7 PC 上需要 16 秒：
n = 400;
m = 1e6;
# Generate data
mat = matrix(runif(m*n),n,m);
# Start timer
tic = proc.time();
# Center each variable
mat = mat - rowMeans(mat);
# Standardize each variable
mat = mat / sqrt(rowSums(mat^2));
# Calculate correlations
cr = tcrossprod(mat);
# Stop timer
toc = proc.time();
# Show the results and the time
show(cr[1:4,1:4]);
show(toc-tic)
上面的 R 代码报告了以下时序：
user system elapsed
31.82 1.98 15.74
我在我的MatrixEQTL 包中使用了这种方法。http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/
有关 R 的各种 BLAS 选项的更多信息，请访问此处：http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/runit.html#large
【讨论】：
不使用任何优化的 BLAS 构建 R，在我的机器 (2.9Ghz i7) 上大约需要 2 分钟。我会用优化的 BLAS 安装 R 并告诉你。 是的，@user1971988，我很好奇这段代码在 BLAS 中的性能。 另外，如果您愿意，本网站的习惯是接受答案。 在使用优化的 BLAS 从源代码重新安装 R 后，我正在尝试复制您的时间。给我几天时间，我会更新我的结果并接受你的回答。 它使用什么方法？【参考方案2】：一些事情。
首先，如果您使用 .C 接口进行外部调用，则默认情况下它会复制所有参数。这就是对象 corr 没有被修改的原因。如果您想避免这种情况，则必须在 .C 调用中设置 DUP=false。但是，通常使用 .C 来修改现有的 R 对象并不是做事的首选方式。相反，您可能想要创建一个新数组并允许外部调用来填充它，就像这样。
corr<-.C("rCorrelationWrapper2", as.double(X), as.integer(dim(X)),
as.double(mu), as.double(sd),
as.integer(Rows), as.integer(Cols),
result=double(p*q))$result
corr<-array(corr,c(p,q))
其次，就编写快速相关函数而言，您应该尝试的第一件事是使用高效的 BLAS 实现来编译 R。这不仅会使您的相关函数更快，还会使您的所有线性代数更快。好的免费候选人是 AMD 的 ACML 或 ATLAS。其中任何一个都能够非常快速地计算相关矩阵。加速不仅仅是并行化——这些库在缓存使用方面也很聪明，并且在程序集级别进行了优化，因此即使只有一个内核，您也会看到很大的改进。 http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acml/ http://math-atlas.sourceforge.net/
最后，如果您真的想编写自己的 C 代码，我建议您使用 openMP 在不同线程之间自动拆分计算，而不是手动完成。但是，对于像矩阵乘法这样基本的东西，最好使用可用的优化库。
【讨论】：
以上是关于使用 C 和并行化在 R 中快速关联的主要内容，如果未能解决你的问题，请参考以下文章
为啥在更多 CPU/内核上的并行化在 Python 中的扩展性如此之差？
R：如何在 R 3.2.1 中使用 lattice 并行化多面板绘图？