在 R 中创建矩阵的更快方法

Question

我有以下我一直在做的事情的例子，这在形式上很简单，但我想检查一下我的代码有哪些潜在的替代方案——如果可能的话，为了更快。这是一个例子：

Time1=Sys.time()
v=rep(c("A","B","C","D","E","F","G","H","I","J","K","L","M","N","O","P"),
     each=1000)
m=matrix(0,ncol=length(v),nrow=length(v))
for (j in 1:length(v)) {
  for(i in 1:length(v)) {
    if (v[j]==v[i]) {
      m[j,i]=1
    } else {
      next
    }
  }
}
Time2=Sys.time()
Time2-Time1
# Time difference of 1.405404 mins

我正在创建一个简单的关系矩阵——其中向量v1 可以被解释为被放置为行和列以及结果相等的矩阵映射。如果它们相等，我们得到m[j,i]=1；如果不相等，m[j,i]=0。正如我所说，我想让这段代码运行得更快。我试图想办法将它编码为apply 函数，但我现在还没有想出来。不过，我想知道除了我所说的之外是否还有其他选择。

编辑：我对文本做了一些更正，并试图澄清问题。

Answer 1

当然。假设您的样本数据不代表真实数据，那么它的工作速度大约快 6 倍：

m2 <- +(outer(v, v, `==`))
all.equal(m, m2)
# [1] TRUE

但是，如果您的真实数据有大量重复，那么 @Sathish 在比较之前删除重复并通过矩阵传播的方法可能会快得多。

Answer 2

我认为@r2evans 的outer 方法是构造矩阵的最简单方法。下面是另一个基本 R 选项，使用 expand.grid

m2 <- matrix(+do.call("==",expand.grid(v,v)),length(v))

Answer 3

我将与来自data.table 的CJ 一起参加比赛。

libary(data.table)
m2 <- matrix(+(CJ(v1 = v,v2 = v,sorted=FALSE)[,ans := v1==v2][,ans]),length(v))
all.equal(m,m2)
#[1] TRUE

Answer 4

如果你有很多零，那么你可以使用稀疏矩阵。您填写匹配的位置，其余位置为零。这会消耗更少的内存，但您只能将其与某些功能一起使用。例如在 glmnet lasso 中使用它。

我认为@r2evans 解决方案是最简洁的，适用于大多数情况。

下面我有一些来自答案的代码，它们中的大多数确实比 OP 的要快

library(microbenchmark)
library(Matrix)
library(data.table)
setDTthreads(threads =1)

f_sw = function(v){
N = length(v)
i = lapply(v,function(i)which(v==i))
j = rep(1:N,times=sapply(i,length))
as.matrix(sparseMatrix(i=unlist(i),j=j,dims=list(N,N)))
}

f_r2evans = function(v){
  m2 <- +(outer(v, v, `==`))
  return(m2)
}

f_IanCampbell = function(v){ 
  matrix(+(CJ(v1 = v,v2 = v,sorted=FALSE)[,ans := v1==v2][,ans]),length(v))
}

microbenchmark(f_IanCampbell(v),f_sw(v),f_r2evans(v),times=5)

Unit: seconds
             expr       min        lq      mean    median        uq      max
 f_IanCampbell(v) 10.820616 11.325422 12.544146 12.983926 13.126655 14.46411
          f_sw(v)  7.014364  7.228585  8.206858  7.745741  8.877425 10.16818
     f_r2evans(v)  9.117405  9.519443  9.996789  9.896823 10.288586 11.16169
 neval cld
     5   b
     5  a 
     5  a

Answer 5

可能是这样的：

我们可以在向量上执行此操作，然后逐行逐列地将矩阵复制到 1000 次，而不是检查所有可能值的相等性。它将给出相同的结果。此代码不维护列和行顺序。但是，使用行名和列名，我们可以验证答案是否正确。

我使用t()，因为列绑定比行绑定快。

system.time({
  v1 <- c("A","B","C","D","E","F","G","H","I","J","K","L","M","N","O","P")
  m1 <- sapply(v1, function(x) as.integer(v1 == x))
  rownames(m1) <- colnames(m1)
  m1 <- do.call('cbind', mget(rep('m1', 1000)))
  m1 <- t(m1)
  m1 <- do.call('cbind', mget(rep('m1', 1000)))
  m1 <- t(m1)
}) 
# user  system elapsed 
# 9.32    0.50    9.84 


dim(m1)
# [1] 16000 16000

---

另一种方法：

这个（硬编码）不会进行任何比较，但我们会根据可能发生的情况来创建值，方法是将向量与其自身的值进行比较。

可以通过使用eval(parse(text=paste())) 构造来改进。

system.time({
  m4 <- matrix(data = 
           c(
             c(rep(1, 1000), rep(0, 15000)),
             c(rep(0, 1000), rep(1, 1000), rep(0, 14000)),
             c(rep(0, 2000), rep(1, 1000), rep(0, 13000)),
             c(rep(0, 3000), rep(1, 1000), rep(0, 12000)),
             c(rep(0, 4000), rep(1, 1000), rep(0, 11000)),
             c(rep(0, 5000), rep(1, 1000), rep(0, 10000)),
             c(rep(0, 6000), rep(1, 1000), rep(0, 9000)),
             c(rep(0, 7000), rep(1, 1000), rep(0, 8000)),
             c(rep(0, 8000), rep(1, 1000), rep(0, 7000)),
             c(rep(0, 9000), rep(1, 1000), rep(0, 6000)),
             c(rep(0, 10000), rep(1, 1000), rep(0, 5000)),
             c(rep(0, 11000), rep(1, 1000), rep(0, 4000)),
             c(rep(0, 12000), rep(1, 1000), rep(0, 3000)),
             c(rep(0, 13000), rep(1, 1000), rep(0, 2000)),
             c(rep(0, 14000), rep(1, 1000), rep(0, 1000)),
             c(rep(0, 15000), rep(1, 1000))), nrow = 16000, ncol = 16000)
})
# user  system elapsed 
# 0.72    0.93    1.82 

注意： 正如@r2evans 所说，如果 OP 的样本数据不代表真实数据，这将不起作用