Запуск cor () (или любого другого варианта) над разреженной матрицей в R

Question

При попытке запустить функцию cor() на разреженных матрицах (типа dgCMatrix или dgTMatrix) я получаю следующую ошибку:

Error in cor(x) : supply both 'x' and 'y' or a matrix-like 'x'

Преобразование моей матрицы в плотную будет очень неэффективным. Есть ли простой способ рассчитать эту корреляцию (без цикла всех пар?).

Спасибо

• Рон

Answer 1

EDITED ANSWER - оптимизирован для использования памяти и скорости.

Ваша ошибка логична, поскольку разреженная матрица не распознается функцией cor как матрица, и в пакете Matrix отсутствует метод -yet- no для корреляций.

Я не знаю ни одной функции, которая позволила бы вам рассчитать это, но вы можете легко вычислить это самостоятельно, используя матричные операторы, доступные в пакете Matrix:

sparse.cor <- function(x){
  n <- nrow(x)
  m <- ncol(x)
  ii <- unique(x@i)+1 # rows with a non-zero element

  Ex <- colMeans(x)
  nozero <- as.vector(x[ii,]) - rep(Ex,each=length(ii))        # colmeans

  covmat <- ( crossprod(matrix(nozero,ncol=m)) +
              crossprod(t(Ex))*(n-length(ii))
            )/(n-1)
  sdvec <- sqrt(diag(covmat))
  covmat/crossprod(t(sdvec))
}

covmat - это ваша матрица дисперсии-ковариации, так что вы можете вычислить и ее. Расчет основан на выборе строк, где хотя бы один элемент не равен нулю. к кросс-произведению этого вы добавляете числа Колм, умноженные на количество строк с нулем. Это эквивалентно

(X - E [X]) раз (X - E [X]) транспонировано

Разделите на n-1, и у вас будет матрица дисперсии-ковариации. Остальное легко.

Контрольный пример:

X <- sample(0:10,1e8,replace=T,p=c(0.99,rep(0.001,10)))
xx <- Matrix(X,ncol=5)

> system.time(out1 <- sparse.cor(xx))
   user  system elapsed 
   0.50    0.09    0.59 
> system.time(out2 <- cor(as.matrix(xx)))
   user  system elapsed 
   1.75    0.28    2.05 
> all.equal(out1,out2)
[1] TRUE

Answer 2

Это то, что я в конечном итоге использовал. Спасибо @Joris за помощь!

Моя x матрица довольно большая. Предполагая, что это размер n * p, n=200k и p=10k в моем случае.

Хитрость заключается в том, чтобы поддерживать разреженность операций и выполнять вычисления на p * p матрицах вместо p*n x n*p.

Версия 1 более проста, но менее эффективна по времени и памяти, поскольку операция с внешним продуктом стоит дорого:

sparse.cor2 <- function(x){
    n <- nrow(x)

    covmat <- (crossprod(x)-2*(colMeans(x) %o% colSums(x))
              +n*colMeans(x)%o%colMeans(x))/(n-1)

    sdvec <- sqrt(diag(covmat)) # standard deviations of columns
    covmat/crossprod(t(sdvec)) # correlation matrix
}

Версия 2 более эффективна по времени (экономит несколько операций) и по памяти. Все еще требует огромного количества памяти для матрицы p=10k:

sparse.cor3 <- function(x){
    memory.limit(size=10000)
    n <- nrow(x)

    cMeans <- colMeans(x)
    cSums <- colSums(x)

    # Calculate the population covariance matrix.
    # There's no need to divide by (n-1) as the std. dev is also calculated the same way.
    # The code is optimized to minize use of memory and expensive operations
    covmat <- tcrossprod(cMeans, (-2*cSums+n*cMeans))
    crossp <- as.matrix(crossprod(x))
    covmat <- covmat+crossp

    sdvec <- sqrt(diag(covmat)) # standard deviations of columns
    covmat/crossprod(t(sdvec)) # correlation matrix
}

Сравнение времени (sparse.cor - последняя версия @Joris):

> X <- sample(0:10,1e7,replace=T,p=c(0.9,rep(0.01,10)))
> x <- Matrix(X,ncol=10)
> 
> object.size(x)
11999472 bytes
> 
> system.time(corx <- sparse.cor(x))
   user  system elapsed 
   0.50    0.06    0.56 
> system.time(corx2 <- sparse.cor2(x))
   user  system elapsed 
   0.17    0.00    0.17 
> system.time(corx3 <- sparse.cor3(x))
   user  system elapsed 
   0.13    0.00    0.12 
> system.time(correg <-cor(as.matrix(x)))
   user  system elapsed 
   0.25    0.03    0.29 
> all.equal(c(as.matrix(corx)),c(as.matrix(correg)))
[1] TRUE
> all.equal(c(as.matrix(corx2)),c(as.matrix(correg)))
[1] TRUE
> all.equal(c(as.matrix(corx3)),c(as.matrix(correg)))
[1] TRUE

намного больше x матрица:

> X <- sample(0:10,1e8,replace=T,p=c(0.9,rep(0.01,10)))
> x <- Matrix(X,ncol=10)
> object.size(x)
120005688 bytes
> system.time(corx2 <- sparse.cor2(x))
   user  system elapsed 
   1.47    0.07    1.53 
> system.time(corx3 <- sparse.cor3(x))
   user  system elapsed 
   1.18    0.09    1.29 
> system.time(corx <- sparse.cor(x))
   user  system elapsed 
   5.43    1.26    6.71

Answer 3

Ответ был элегантно решен @Ron, но небольшая модификация решения немного чище и также возвращает образец ковариационной матрицы.

sparse.cor4 <- function(x){
    n <- nrow(x)
    cMeans <- colMeans(x)
    covmat <- (as.matrix(crossprod(x)) - n*tcrossprod(cMeans))/(n-1)
    sdvec <- sqrt(diag(covmat)) 
    cormat <- covmat/tcrossprod(sdvec)
    list(cov=covmat,cor=cormat)
}

Упрощение проистекает из этого: с матрицей X x p и матрицей M n x p столбца означает X:

cov(X) = E[(X-M)'(X-M)] = E[X'X - M'X - X'M + M'M] 

M'X = X'M = M'M, which have (i,j) elements = sum(column i) * sum(column j) / n

= n * mean(column i) * mean(column j)

или записано с вектором строки m столбца означает,

= n * m'm

Тогда cov(X) = E[X'X - n m'm]

---

и теперь он чуть быстрее.

> X <- sample(0:10,1e7,replace=T,p=c(0.9,rep(0.01,10)))
> x <- Matrix(X,ncol=10)
> system.time(corx <- sparse.cor(x))
   user  system elapsed 
  1.139   0.196   1.334 
> system.time(corx3 <- sparse.cor3(x))
   user  system elapsed 
  0.194   0.007   0.201 
> system.time(corx4 <- sparse.cor4(x))
   user  system elapsed 
  0.187   0.007   0.194 
> system.time(correg <-cor(as.matrix(x)))
   user  system elapsed 
  0.341   0.067   0.407 
> system.time(covreg <- cov(as.matrix(x)))
   user  system elapsed 
  0.314   0.016   0.330 
> all.equal(c(as.matrix(corx)),c(as.matrix(correg)))
[1] TRUE
> all.equal(c(as.matrix(corx3)),c(as.matrix(correg)))
[1] TRUE
> all.equal(c(as.matrix(corx4$cor)),c(as.matrix(correg)))
[1] TRUE
> all.equal(c(as.matrix(corx4$cov)),c(as.matrix(covreg)))
[1] TRUE