r - 双循环的更快计算?
问题描述
我有一段工作代码需要花费太多时间(几天?)来计算。我有一个 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 的矩阵,其中每个条目是这些操作的结果。这个循环有效,但有什么想法可以提高效率吗?谢谢
解决方案
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 元素的行有意义
==> 我们可以跳过所有m[, j]
只包含 0 的计算,所以:
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)
推荐阅读
- git - 如何为 repo 更新全局 git 钩子?
- types - RxJS ids 从类型字符串到类型数字
- php - 在 PHP 中将 XML 对象转换为数组不起作用
- sql - Postgresql,查询以计算特定性别的机器使用次数
- python - ProcessPoolExecutor pass multiple arguments
- c++ - 为什么这些范围之外的字符也会被输出文件中的加法/减法改变?
- python - 需要从与用户表相关的日期表中获取值
- visual-studio - VS2019安装SQL Server组件查看网页时出现err_connection_reset
- docker - GKE - 在运行时绕过 Pod LoadBalancer(Pod 的外部 IP)到 Pod 的容器的 IP 以用于 WebSocket 目的
- epub - 在 epubcheck 中验证 epub 时出错,它说这个