R：向量化循环内对象修改

Question

我想知道是否可以向量化当前使用循环的函数。

给定一个示例矩阵：

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() 的性能对比改进的循环函数（详见下文）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

Answer 1

这是另一种存储行和列总和并在选择行和列后更新的选项：

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：向量化循环内对象修改

R: Vectorize in-loop object modification

performance

for-loop

r

matrix

vectorization