0

我正在做一些需要我重复计算大型方阵的元素的事情。该过程涉及读取存储在另一个矩阵中的数据，然后计算矩阵元素。目前我正在使用双for循环来执行此操作。

library(matrixcalc)

data <- matrix(nrow=3,ncol=1000)

for(x in 1:ncol(data)){
   for(y in 1:ncol(data)){
       matrix[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2))
   }
}

问题是这非常慢，因为我的矩阵非常大。这个过程最快的替代方法是什么？

4

4 回答 4

3

短且非常快：

mat <- exp(-as.matrix(dist(t(data))))

我还建议将该fields::rdist函数作为计算欧几里德距离矩阵的更快替代方法dist，因此如果加载包不是问题，请考虑：

library(fields)
mat <- exp(-rdist(t(data)))

为了让您了解速度改进：

data <- matrix(runif(3000), nrow=3, ncol=1000)

OP <- function(data) {
  require(matrixcalc)
  mat <- matrix(0, ncol(data), ncol(data))
  for(x in 1:ncol(data)){
    for(y in 1:ncol(data)){
      mat[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2))
    }
  }
  mat
}

flodel1 <- function(data) exp(-as.matrix(dist(t(data))))
flodel2 <- function(data) {
  require(fields)
  exp(-rdist(t(data)))
}

system.time(res1 <- OP(data))
#   user  system elapsed 
# 22.708   2.080  24.602 
system.time(res2 <- flodel1(data))
#   user  system elapsed 
#  0.112   0.025   0.136 
system.time(res3 <- flodel2(data))
#   user  system elapsed 
#  0.048   0.000   0.049

（请注意，在 和 的情况下OP，flodel2这些运行时不包括在测试之前加载的包的加载。）

于 2013-07-04T11:44:11.880 回答
2

这应该快得多：

nc <- ncol(data)

mat <- diag(nc)

for(x in 2:nc){
   for(y in 1:x){
       mat[x, y] <- exp(-(sum((data[ , x] - data[ , y])^2) ^ .5))
   }
}

mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
于 2013-07-04T11:03:15.087 回答
1

R 语言使用列主序数组。更改 for 循环顺序可以提高性能。因为这样，您可以以更连续的形式访问内存，从而实现 cpu-cache 优势。

 for(y in 1:dim) //outer is y now
 {
    for(x in 1:dim) //now x is count inside
    {
        matrix[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2))
    }
 }

你的“矩阵”是二维数组对吗？

如果您需要更快的速度，您可以展开一些内部循环以减少 cpu 的分支负载和更好的缓存/预取。

 for(y in 1:dim) 
 {
    for(x in 1:(dim/8)) //lets imagine dimension is a multiple of 8
    {
        matrix[x,y]=exp(-entrywise.norm(data[,x]-data[,y],2))
        matrix[x+1,y]=exp(-entrywise.norm(data[,x+1]-data[,y],2))
        matrix[x+2,y]=exp(-entrywise.norm(data[,x+2]-data[,y],2))
        matrix[x+3,y]=exp(-entrywise.norm(data[,x+3]-data[,y],2))
        matrix[x+4,y]=exp(-entrywise.norm(data[,x+4]-data[,y],2))
        matrix[x+5,y]=exp(-entrywise.norm(data[,x+5]-data[,y],2))
        matrix[x+6,y]=exp(-entrywise.norm(data[,x+6]-data[,y],2))
        matrix[x+7,y]=exp(-entrywise.norm(data[,x+7]-data[,y],2))
    }
 }
于 2013-07-04T11:00:07.270 回答
1

您可以使用colSums而不是内部循环。根据@Sven Hohenstein 的回答：

nc <- ncol(data)

mat <- diag(nc)

for(x in 2:nc){
  mat[x, 1:x] <- exp(-(colSums((data[ , 1:x] - data[ ,x])^2) ^ .5))
}

mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
于 2013-07-04T11:35:17.523 回答