双循环的更快计算？

Question

我有一段工作代码需要花费太多小时（几天？）来计算。 我有一个 1 和 0 的稀疏矩阵，我需要以所有可能的组合从任何其他行中减去每一行，将结果向量乘以另一个向量，最后平均其中的值以获得我需要的单个标量插入矩阵。我所拥有的是：

m <- matrix( 
c(0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0), nrow=4,ncol=4,
byrow = TRUE)   

b <- c(1,2,3,4)

for (j in 1:dim(m)[1]){
 for (i in 1:dim(m)[1]){
    a <- m[j,] - m[i,]
    a[i] <- 0L
    a[a < 0] <- 0L
    c <- a*b
    d[i,j] <- mean(c[c > 0])
 }
}

所需的输出是具有相同维度 m 的矩阵，其中每个条目是这些操作的结果。 这个循环有效，但有什么想法可以提高效率吗？谢谢

Answer 1

我愚蠢的解决方案是使用apply 或sapply 函数，而不是for 循环来进行迭代：

sapply(1:dim(m)[1], function(k) {z <- t(apply(m, 1, function(x) m[k,]-x)); diag(z) <- 0; z[z<0] <- 0; apply(t(apply(z, 1, function(x) x*b)),1,function(x) mean(x[x>0]))})

我试图比较你的解决方案，这在我的计算机上运行时间方面，你的需要

t1 <- Sys.time()
d1 <- m
for (j in 1:dim(m)[1]){
  for (i in 1:dim(m)[1]){
    a <- m[j,] - m[i,]
    a[i] <- 0L
    a[a < 0] <- 0L
    c <- a*b
    d1[i,j] <- mean(c[c > 0])
  }
}
Sys.time()-t1

您需要Time difference of 0.02799988 secs。对我来说，它减少了一点但不会太多，即Time difference of 0.01899815 secs，当你运行时

t2 <- Sys.time()
d2 <- sapply(1:dim(m)[1], function(k) {z <- t(apply(m, 1, function(x) m[k,]-x)); diag(z) <- 0; z[z<0] <- 0; apply(t(apply(z, 1, function(x) x*b)),1,function(x) mean(x[x>0]))})
Sys.time()-t2

你可以在自己的电脑上试一试，矩阵更大，祝你好运！

Answer 2

1) 创建测试稀疏矩阵：

nc <- nr <- 100
p <- 0.001
require(Matrix)
M <- Matrix(0L, nr, nc, sparse = T) # 0 matrix
n1 <- ceiling(p * (prod(dim(M)))) # 1 count
M[1:n1] <- 1L # fill only first column, to approximate max non 0 row count
# (each row has at maximum 1 positive element)
sum(M)/(prod(dim(M)))

b <- 1:ncol(M)

sum(rowSums(M))

所以，如果给定的比例是正确的，那么我们最多有 10 行包含非 0 元素

基于这一事实和您提供的计算：

# a <- m[j, ] - m[i, ]
# a[i] <- 0L
# a[a < 0] <- 0L
# c <- a*b
# mean(c[c > 0])

我们可以看到结果只对m[, j] 具有至少 1 个非 0 元素的行有意义

==> 我们可以跳过所有只包含 0 的 m[, j] 的计算，所以：

minem <- function() { # write as function
  t1 <- proc.time() # timing
  require(data.table)
  i <- CJ(1:nr, 1:nr) # generate all combinations
  k <- rowSums(M) > 0L # get index where at least 1 element is greater that 0
  i <- i[data.table(V1 = 1:nr, k), on = 'V1'] # merge
  cat('at moust', i[, sum(k)/.N*100], '% of rows needs to be calculated \n')
  i[k == T, rowN := 1:.N] # add row nr for 0 subset
  i2 <- i[k == T] # subset only those indexes who need calculation
  a <- M[i2[[1]],] - M[i2[[2]],] # operate on all combinations at once
  a <- drop0(a) # clean up 0

  ids <- as.matrix(i2[, .(rowN, V2)]) # ids for 0 subset
  a[ids] <- 0L # your line: a[i] <- 0L
  a <- drop0(a) # clean up 0

  a[a < 0] <- 0L # the same as your line
  a <- drop0(a) # clean up 0

  c <- t(t(a)*b) # multiply each row with vector
  c <- drop0(c) # clean up 0

  c[c < 0L] <- 0L # for mean calculation
  c <- drop0(c) # clean up 0

  r <- rowSums(c)/rowSums(c > 0L) # row means
  i[k == T, result := r] # assign results to data.table
  i[is.na(result), result := NaN] # set rest to NaN
  d2 <- matrix(i$result, nr, nr, byrow = F) # create resulting matrix
  t2 <- proc.time() # timing
  cat(t2[3] - t1[3], 'sec \n')
  d2
}
d2 <- minem()
# at most 10 % of rows needs to be calculated 
# 0.05 sec 

如果结果匹配，则测试较小的示例

d <- matrix(NA, nrow(M), ncol(M))
for (j in 1:dim(M)[1]) {
  for (i in 1:dim(M)[1]) {
    a <- M[j, ] - M[i, ]
    a[i] <- 0L
    a[a < 0] <- 0L
    c <- a*b
    d[i, j] <- mean(c[c > 0])
  }
}
all.equal(d, d2)

我们可以得到你真实数据大小的结果吗？：

# generate data:
nc <- nr <- 6663L
b <- 1:nr
p <- 0.0001074096 # proportion of 1s
M <- Matrix(0L, nr, nc, sparse = T) # 0 matrix
n1 <- ceiling(p * (prod(dim(M)))) # 1 count
M[1:n1] <- 1L

object.size(as.matrix(M))/object.size(M)
# storing this data in usual matrix uses 4000+ times more memory

# calculation:
d2 <- minem()
# at most 71.57437 % of rows needs to be calculated 
# 28.33 sec 

所以你需要将你的矩阵转换为稀疏矩阵

M <- Matrix(m, sparse = T)