Быстрая корреляция в R с использованием C и распараллеливания
Мой проект на сегодня заключался в том, чтобы написать процедуру быстрой корреляции на R с использованием базового набора навыков, который у меня есть. Я должен найти корреляцию между почти 400 переменными, каждая из которых имеет почти миллион наблюдений (то есть матрица размером p=1 мм строк и n=400 столбцов).
Собственная корреляционная функция R занимает почти 2 минуты для 1-миллиметровых строк и 200 наблюдений на переменную. Я не проводил 400 наблюдений за колонку, но думаю, это займет почти 8 минут. У меня есть менее 30 секунд, чтобы закончить это.
Следовательно, я хочу делать что-то.
1 - написать простую корреляционную функцию в C и применить ее в блоках параллельно (см. Ниже).
2 - Блоки - разбить корреляционную матрицу на три блока (верхний левый квадрат размера K*K, нижний правый квадрат размера (pK) (pK) и верхняя правая прямоугольная матрица размера K (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, как мне загрузить общий объект в каждое ядро для параллельных заданий на каждом ядре на этапе инициализации (а не как я это сделал)?
Вопрос № 3: Что означает эта ошибка? Мне нужны некоторые указатели, и я хотел бы отладить это сам.
Вопрос № 4: Существует ли быстрый способ вычисления корреляции по матрицам 1 ММ на 400, менее чем за 30 секунд?
Когда я запускаю МЕТОД 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));
}
2 ответа
Используя быстрый 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/
Более подробная информация о различных опциях BLAS для R доступна здесь:
http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/runit.html
Несколько вещей.
Во-первых, если вы используете интерфейс.C для внешних вызовов, то по умолчанию он создает копии всех аргументов. Вот почему объект corr не изменяется. Если вы хотите избежать этого, вы должны установить DUP=false в вызове.C. Однако в общем случае использование.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))
Во-вторых, что касается написания быстрой корреляционной функции, первое, что вы должны попробовать - это скомпилировать R с эффективной реализацией BLAS. Это не только ускорит вашу корреляционную функцию, но и ускорит всю вашу линейную алгебру. Хорошие бесплатные кандидаты - ACML от AMD или ATLAS. Любой из них сможет очень быстро вычислить матрицы корреляции. Ускорение - это больше, чем просто распараллеливание - эти библиотеки также хорошо разбираются в использовании кеша и оптимизируются на уровне сборки, поэтому даже с одним ядром вы увидите значительное улучшение. http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acml/ http://math-atlas.sourceforge.net/
Наконец, если вы действительно хотите написать свой собственный код на C, я бы предложил использовать openMP для автоматического разделения вычислений между различными потоками, а не делать это вручную. Но для чего-то такого базового, как умножение матриц, вероятно, лучше использовать доступную оптимизированную библиотеку.