使用C和并行化在R中快速相关

我今天的项目是使用我拥有的基本技能在R中编写一个快速关联例程。我必须找到几乎400个变量之间的相关性，每个变量有近百万个观测值（即大小为p = 1MM行且n = 400个cols的矩阵）。

R的原生相关函数对于1MM行需要近2分钟，每个变量需要200个观察值。我没有为每列运行400次观察，但我的猜测是需要将近8分钟。我完成它的时间不到30秒。

因此，我想做的事情。

1 – 在C中编写一个简单的相关函数，并将其平行应用于块中（见下文）。

2 – 块 – 将相关矩阵分成三个块（大小为K * K的左上方，大小为右下方（pK） （pK），以及大小为K （pK）的右上方矩形矩阵 ）。这涵盖了相关矩阵corr中的所有单元，因为我只需要上三角形。

3 – 使用降雪并行地通过.C调用运行Cfunction。

 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，如何在init步骤中为每个核心上的并行作业将共享对象加载到每个核心（而不是我如何完成）？

问题3：这个错误是什么意思？ 我需要一些指针，我很乐意自己调试。

问题4：在不到30秒的时间内，是否有一种快速计算矩阵1MM乘以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  #include  #include  #include  #include  // 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("\np: %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("\ni: %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)); }

使用快速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

一些东西。

首先，如果您使用.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或ATLAS的ACML。其中任何一个都能够非常快速地计算相关矩阵。加速不仅仅是并行化 - 这些库对于缓存使用情况也非常聪明，并且在汇编级别进行了优化，因此即使只有一个核心，您也会看到很大的改进。 http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acml/ http://math-atlas.sourceforge.net/

最后，如果你真的想编写自己的C代码，我建议使用openMP自动在不同的线程之间分配计算，而不是手工完成。但是，对于像矩阵乘法这样基本的东西，使用可用的优化库可能会更好。