将函数应用于R中的距离矩阵

Question

这个问题今天出现在操纵邮件列表中.

http://groups.google.com/group/manipulatr/browse_thread/thread/fbab76945f7cba3f

我在改写.

给定距离矩阵(用其计算dist)将函数应用于距离矩阵的行.

码:

library(plyr)
N <- 100
a <- data.frame(b=1:N,c=runif(N))
d <- dist(a,diag=T,upper=T)
sumd <- adply(as.matrix(d),1,sum)

问题是要按行应用函数,你必须存储整个矩阵(而不仅仅是下三角形部分.因此它对大型矩阵使用了太多内存.在我的计算机中,对于大小为10000的矩阵,它会失败.

有任何想法吗？

Answer 1

首先,对于尚未见过这一点的人,我强烈建议您阅读r-wiki上有关代码优化的文章.

这是另一个没有使用的版本ifelse(这是一个相对较慢的功能):

noeq.2 <- function(i, j, N) {
    i <- i-1
    j <- j-1
    x <- i*(N-1) - (i-1)*((i-1) + 1)/2 + j - i
    x2 <- j*(N-1) - (j-1)*((j-1) + 1)/2 + i - j
    idx <- i < j
    x[!idx] <- x2[!idx]
    x[i==j] <- 0
    x
}

和我的笔记本电脑上的时间:

> N <- 1000
> system.time(sapply(1:N, function(i) sapply(1:N, function(j) noeq(i, j, N))))
   user  system elapsed 
  51.31    0.10   52.06 
> system.time(sapply(1:N, function(j) noeq.1(1:N, j, N)))
   user  system elapsed 
   2.47    0.02    2.67 
> system.time(sapply(1:N, function(j) noeq.2(1:N, j, N)))
   user  system elapsed 
   0.88    0.01    1.12 

lapply比sapply更快:

> system.time(do.call("rbind",lapply(1:N, function(j) noeq.2(1:N, j, N))))
   user  system elapsed 
   0.67    0.00    0.67 

Answer 2

我的解决方案是在给定行和矩阵大小的情况下获取距离向量的索引。我从codeguru那里得到了这个

int Trag_noeq(int row, int col, int N)
{
   //assert(row != col);    //You can add this in if you like
   if (row<col)
      return row*(N-1) - (row-1)*((row-1) + 1)/2 + col - row - 1;
   else if (col<row)
      return col*(N-1) - (col-1)*((col-1) + 1)/2 + row - col - 1;
   else
      return -1;
}

转换为 R 后，假设索引从 1 开始，并假设我得到的是下三矩阵而不是上三矩阵。
编辑：使用 rcs 贡献的矢量化版本

noeq.1 <- function(i, j, N) {
    i <- i-1
    j <- j-1
    ix <- ifelse(i < j,
                 i*(N-1) - (i-1)*((i-1) + 1)/2 + j - i,
                 j*(N-1) - (j-1)*((j-1) + 1)/2 + i - j) * ifelse(i == j, 0, 1)
    ix
}

## To get the indexes of the row, the following one liner works:

getrow <- function(z, N) noeq.1(z, 1:N, N)

## to get the row sums

getsum <- function(d, f=sum) {
    N <- attr(d, "Size")
    sapply(1:N, function(i) {
        if (i%%100==0) print(i)
        f(d[getrow(i,N)])
    })
}

所以，举个例子：

sumd2 <- getsum(d)

对于矢量化之前的小矩阵，这比 as.matrix 慢得多。但矢量化后速度慢了大约 3 倍。在 Intel Core2Duo 2ghz 中，按行应用大小为 10000 的矩阵的总和只需 100 多秒。as.matrix 方法失败。谢谢rcs！