Быстрая корреляция в R с использованием C и распараллеливания

Мой проект на сегодня состоял в том, чтобы написать программу быстрой корреляции на R с использованием имеющегося у меня базового набора навыков. Мне нужно найти корреляцию между почти 400 переменными, каждая из которых имеет почти миллион наблюдений (т.е. матрица размером p = 1MM строк и n = 400 столбцов).

Собственная корреляционная функция R занимает почти 2 минуты для 1MM строк и 200 наблюдений на переменную. Я не проводил 400 наблюдений на столбец, но предполагаю, что это займет почти 8 минут. У меня меньше 30 секунд, чтобы закончить.

Следовательно, я хочу что-то делать.

1 - напишите простую корреляционную функцию на C и примените ее параллельно блоками (см. Ниже).

2 - Блоки - разделите матрицу корреляции на три блока (верхний левый квадрат размера K * K, нижний правый квадрат размера (pK) (pK) и верхний правый прямоугольный квадрат размера K ( рК)). Это покрывает все ячейки в корреляционной матрице 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. Есть ли быстрый способ вычисления корреляции по матрицам 1MM на 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));
}

c r correlation snowfall
person user1971988    schedule 23.09.2013    source источник
comment
@MartinMorgan - собственная функция cor в R (на основе имеющейся у меня сборки) занимает больше времени, как я упоминал выше. Я использую предложение Андрея ниже, и это занимает около 2 минут для 1MM на 400 варов. Буду обновлять.   -  person user1971988    schedule 24.09.2013

Ответы (2)


arrow_upward
13
arrow_downward

Используя быстрый BLAS (через Revolution R или Goto BLAS), вы можете быстро вычислить все эти корреляции в R без написания кода C. На моем ПК Intel i7 первого поколения это занимает 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#large

person Andrey Shabalin    schedule 23.09.2013
comment
Без сборки R с использованием какого-либо оптимизированного BLAS на моем компьютере (2,9 ГГц i7) это займет около 2 минут. Я установлю R с оптимизированным BLAS и дам вам знать. - person user1971988; 24.09.2013
comment
Да, @ user1971988, мне было бы интересно узнать, насколько эффективен этот код для вас с BLAS. - person Andrey Shabalin; 25.09.2013
comment
Кроме того, на этом сайте принято принимать ответ, если он вам нравится. - person Andrey Shabalin; 25.09.2013
comment
Я пытаюсь воспроизвести ваше время после переустановки R из исходного кода с использованием оптимизированного BLAS. Дайте мне пару дней, и я обновлю свои результаты и приму ваш ответ. - person user1971988; 25.09.2013
comment
Какой метод он использует? - person user1436187; 17.09.2014
comment
Разве у вас не должно быть cr = tcrossprod (mat) / nrow (mat) (альтернативно cr = tcrossprod (mat) / (nrow (mat) - 1) в зависимости от оценки)? - person user2280549; 30.11.2018
comment
@ user2280549, нет, не надо. Код правильный, как написано, и вы можете проверить его, чтобы подтвердить. - person Andrey Shabalin; 02.12.2018
comment
@AndreyShabalin Я сделал. n = 1000; m = 10; set.seed (12345); mat = матрица (runif (m * n), n, m); corr1 ‹- кор (мат) [1: 4, 1: 4]; mat = mat - rowMeans (мат); мат = мат / sqrt (rowSums (mat ^ 2)); corr2 = tcrossprod (мат); all (corr1 == corr2) Вы можете проверить самостоятельно, что это не соответствует стандартной реализации. Даже размеры вашей матрицы неверны. - person user2280549; 02.12.2018
comment
Я сопоставляю строки, а не столбцы. - person Andrey Shabalin; 03.12.2018

arrow_upward
4
arrow_downward

Несколько вещей.

Во-первых, если вы используете интерфейс .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 для автоматического разделения вычислений между разными потоками, а не делать это вручную. Но для чего-то столь же простого, как умножение матриц, вероятно, лучше использовать доступную оптимизированную библиотеку.

person mrip    schedule 23.09.2013