R: Vectorize in-loop object modification

user_11078618 0
Viewed 89

I want to know if it is possible to vectorize a function that currently uses a loop.

Given an example matrix:

m <- matrix(c(0,2,1,0,0,2,2,1,0), nrow = 3)
row.names(m) <- colnames(m) <- c("apple", "orange", "pear")

I want to find the item with the minimum value for the ratio of rowSums() to rowSums() + colSums(). Whichever item is identified as the minimum is then appended to a vector z and removed from m, and the process repeats until all items have been ordered in z.

The following loop works fine:

    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)
    }

But this loop is slow, given a large enough problem.

I want to know if it is possible to create a vectorized equivalent to this loop function. If this isn't possible, some ideas to improve speed are welcome.

Important

It is important to understand that the removal of items from m will affect subsequent ratio values. For example, the initial ratio values for m are:

apple orange   pear 
   0.4    0.6    0.5

In this case, in the first iteration, apple would be removed from m and appended to z.

In the next iteration, the ratio values for the remaining items are:

   orange      pear 
0.3333333 0.6666667

So you can see that the ratio value is dependent on the items remaining in tmp.mat.

Update

Performance of loop.function() vs improved loop function (detailed below) lf2() vs Rcpp function 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
vectorization r for-loop matrix performance
1 Answers

Here is another option to store the row and column sums and update once a row and column has been chosen:

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]
}

check:

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

timing on actual dimensions and iterations:

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

More speed gains by coding this in 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;
}
')

timings using Rcpp using a 15 x 15 matrix:

system.time(replicate(1e6, recmin(m)))
#   user  system elapsed 
#   6.17    0.02    6.20
user_19894800
Related