我想知道是否可以向量化当前使用循环的函数。
给定一个示例矩阵:
m <- matrix(c(0,2,1,0,0,2,2,1,0), nrow = 3)
row.names(m) <- colnames(m) <- c("apple", "orange", "pear")
我想找到 rowSums()
的比例值最小的项目至 rowSums() + colSums()
.然后将被标识为最小值的项目附加到向量 z
并从 m
中删除,然后重复该过程,直到在 z
中订购了所有商品为止.
下面的循环工作正常:
loop.function <- function(mat){
nt <- nrow(mat)
z <- rep(NA, nt)
tmp.mat <- mat
for (i in 1:(nt - 1)) {
## ratio value
rv <- rowSums(tmp.mat) / (rowSums(tmp.mat) + colSums(tmp.mat))
## minimum of the ratio values (edited following comment)
min.rv <- which.min(rv)
## append item with minimum ratio value to ith position of z
z[i] <- names(rv)[min.rv]
## remove item appended to z from matrix
tmp.mat <- tmp.mat[-min.rv,-min.rv, drop = FALSE]
}
## append last remaining item of matrix to last position of z
z[nt] <- row.names(tmp.mat)
return(z)
}
但是如果问题足够大,这个循环会很慢。
我想知道是否可以创建一个与此循环函数等效的向量化函数。如果这不可能,欢迎提出一些提高速度的想法。
重要
重要的是要了解从 m
中删除项目会影响后续的比率值。例如,m
的初始比率值是:
apple orange pear
0.4 0.6 0.5
在这种情况下,在第一次迭代中,apple
将从 m
中删除并附加到 z
.
在下一次迭代中,剩余项目的比率值为:
orange pear
0.3333333 0.6666667
因此您可以看到比率值取决于 tmp.mat
中剩余的项目.
更新
loop.function()
的表现vs改进的loop函数(详见下文)lf2()
对比Rcpp
函数 recmin()
:
Unit: microseconds
expr min lq mean median uq max neval cld
loop.function(m) 32.801 33.601 36.33707 34.201 34.9510 75.601 100 c
lf2(m) 20.800 21.701 24.81191 22.151 22.6505 82.200 100 b
recmin(m) 1.601 2.102 2.85100 2.701 3.1000 20.301 100 a
最佳答案
这是另一种存储行和列总和并在选择行和列后更新的选项:
lf2 <- function(m) {
nr <- nrow(m)
res <- integer(nr)
rs <- rowSums(m)
cs <- colSums(m)
for (i in 1L:(nr - 1L)) {
mrv <- which.max(cs / rs)
res[i] <- mrv
rs <- rs - m[, mrv]
cs <- cs - m[mrv,]
cs[mrv] <- -Inf
rs[mrv] <- Inf
}
res[nr] <- which(cs!=-Inf)
rownames(m)[res]
}
检查:
m <- matrix(c(0,2,1,0,0,2,2,1,0), nrow = 3)
row.names(m) <- colnames(m) <- c("apple", "orange", "pear")
identical(loop.function(m), lf2(m))
#[1] TRUE
system.time(replicate(1e5, loop.function(m)))
# user system elapsed
# 3.49 0.00 3.50
system.time(replicate(1e5, lf2(m)))
# user system elapsed
# 1.75 0.00 1.75
实际维度和迭代的时间安排:
set.seed(0L)
n <- 15L
m <- matrix(sample(0L:2L, n*n, TRUE), nrow=n)
rownames(m) <- colnames(m) <- 1L:n
#system.time(replicate(1e5, loop.function(m)))
#Error in z[i] <- names(rv)[min.rv] : replacement has length zero
system.time(replicate(1e5, lf2(m)))
# user system elapsed
# 6.16 0.00 6.16
system.time(replicate(1e6, lf2(m)))
# user system elapsed
# 71.35 0.17 71.55
通过在 Rcpp
中编码可以提高速度:
library(Rcpp)
cppFunction('
IntegerVector recmin(NumericMatrix m) {
int n = m.nrow(), i, j, mrv;
NumericVector rs(n), cs(n);
IntegerVector res(n);
for (i=0; i<n; i++) {
rs[i] = 0.0;
for (j=0; j<n; j++) {
rs[i] += m(i,j);
}
}
for (j=0; j<n; j++) {
cs[j] = 0.0;
for (i=0; i<n; i++) {
cs[j] += m(i,j);
}
}
for (i=0; i<n; i++) {
mrv = n;
for (j=0; j<n; j++) {
if (cs[j] != R_NegInf) {
if (mrv == n) {
mrv = j;
} else if (cs[j] / rs[j] > cs[mrv] / rs[mrv]) {
mrv = j;
}
}
}
res[i] = mrv + 1;
for (j=0; j<n; j++) {
rs[j] -= m(j, mrv);
}
for (j=0; j<n; j++) {
cs[j] -= m(mrv, j);
}
cs[mrv] = R_NegInf;
}
return res;
}
')
使用 Rcpp
使用 15 x 15 矩阵的计时:
system.time(replicate(1e6, recmin(m)))
# user system elapsed
# 6.17 0.02 6.20
关于R:向量化循环内对象修改,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/62759423/