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