Fastest way to find the index of the maximum of each row

Question

I am trying to find an optimal way for finding the index of the maximum value in each row. The problem is that I cannot find a really efficient way in doing it. An example:

Dummy <- matrix(runif(500000000,0,3), ncol = 10000)
> system.time(max.col(Dummy, "first"))
   user  system elapsed 
  5.532   0.075   5.599 
> system.time(apply(Dummy,1,which.max))
   user  system elapsed 
 14.638   0.210  14.828 
> system.time(rowRanges(Dummy))
   user  system elapsed 
  2.083   0.029   2.109 

My main question is, why is it more than 2 times so slow to calculate the indices of the max value in comparison with calculating the max and the min with the rowRanges function. Is there a way how I can improve the performance of calculating the index of the max of each row?

Answer 1

Expanding on krlmlr's answer, some benchmarks:

On dataset:

set.seed(007); Dummy <- matrix(runif(50000000,0,3), ncol = 1000)

maxCol_R is an R by-column loop, maxCol_col is a C by-column loop, maxCol_row is a C by-row loop.

microbenchmark::microbenchmark(max.col(Dummy, "first"), maxCol_R(Dummy), maxCol_col(Dummy), maxCol_row(Dummy), times = 30)
#Unit: milliseconds
#                    expr        min         lq     median         uq       max neval
# max.col(Dummy, "first") 1209.28408 1245.24872 1268.34146 1291.26612 1504.0072    30
#         maxCol_R(Dummy) 1060.99994 1084.80260 1099.41400 1154.11213 1436.2136    30
#       maxCol_col(Dummy)   86.52765   87.22713   89.00142   93.29838  122.2456    30
#       maxCol_row(Dummy)  577.51613  583.96600  598.76010  616.88250  671.9191    30
all.equal(max.col(Dummy, "first"), maxCol_R(Dummy))
#[1] TRUE
all.equal(max.col(Dummy, "first"), maxCol_col(Dummy))
#[1] TRUE
all.equal(max.col(Dummy, "first"), maxCol_row(Dummy))
#[1] TRUE

And the functions:

maxCol_R = function(x)
{
    ans = rep_len(1L, nrow(x))
    mx = x[, 1L]

    for(j in 2:ncol(x)) {
        tmp = x[, j]
        wh = which(tmp > mx)

        ans[wh] = j
        mx[wh] = tmp[wh]
    }

    ans
} 

maxCol_col = inline::cfunction(sig = c(x = "matrix"), body = '
    int nr = INTEGER(getAttrib(x, R_DimSymbol))[0], nc = INTEGER(getAttrib(x, R_DimSymbol))[1]; 
    double *px = REAL(x), *buf = (double *) R_alloc(nr, sizeof(double));
    for(int i = 0; i < nr; i++) buf[i] = R_NegInf;

    SEXP ans = PROTECT(allocVector(INTSXP, nr));
    int *pans = INTEGER(ans);

    for(int j = 0; j < nc; j++) {
        for(int i = 0; i < nr; i++) {
            if(px[i + j*nr] > buf[i]) {
                buf[i] = px[i + j*nr];
                pans[i] = j + 1;
            }
        }
    }

    UNPROTECT(1);
    return(ans);
', language = "C")

maxCol_row = inline::cfunction(sig = c(x = "matrix"), body = '
    int nr = INTEGER(getAttrib(x, R_DimSymbol))[0], nc = INTEGER(getAttrib(x, R_DimSymbol))[1]; 
    double *px = REAL(x), *buf = (double *) R_alloc(nr, sizeof(double));
    for(int i = 0; i < nr; i++) buf[i] = R_NegInf;

    SEXP ans = PROTECT(allocVector(INTSXP, nr));
    int *pans = INTEGER(ans);

    for(int i = 0; i < nr; i++) {
        for(int j = 0; j < nc; j++) {
            if(px[i + j*nr] > buf[i]) {
                buf[i] = px[i + j*nr];
                pans[i] = j + 1;
            }
        }
    }

    UNPROTECT(1);
    return(ans);
', language = "C")

EDIT Jun 10 '16

With slight changes to find the indices of both max and min:

rangeCol = inline::cfunction(sig = c(x = "matrix"), body = '
    int nr = INTEGER(getAttrib(x, R_DimSymbol))[0], nc = INTEGER(getAttrib(x, R_DimSymbol))[1]; 
    double *px = REAL(x), 
           *maxbuf = (double *) R_alloc(nr, sizeof(double)),
           *minbuf = (double *) R_alloc(nr, sizeof(double));
    memcpy(maxbuf, &(px[0 + 0*nr]), nr * sizeof(double));
    memcpy(minbuf, &(px[0 + 0*nr]), nr * sizeof(double));

    SEXP ans = PROTECT(allocMatrix(INTSXP, nr, 2));
    int *pans = INTEGER(ans); 
    for(int i = 0; i < LENGTH(ans); i++) pans[i] = 1;

    for(int j = 1; j < nc; j++) {
        for(int i = 0; i < nr; i++) {
            if(px[i + j*nr] > maxbuf[i]) {
                maxbuf[i] = px[i + j*nr];
                pans[i] = j + 1;
            }
            if(px[i + j*nr] < minbuf[i]) {
                minbuf[i] = px[i + j*nr];
                pans[i + nr] = j + 1;
            }
        }
    }

    UNPROTECT(1);
    return(ans);
', language = "C")

set.seed(007); m = matrix(sample(24) + 0, 6, 4)
m
#     [,1] [,2] [,3] [,4]
#[1,]   24    7   23    6
#[2,]   10   17   21   11
#[3,]    3   22   20   14
#[4,]    2   18    1   15
#[5,]    5   19   12    8
#[6,]   16    4    9   13
rangeCol(m)
#     [,1] [,2]
#[1,]    1    4
#[2,]    3    1
#[3,]    2    1
#[4,]    2    3
#[5,]    2    1
#[6,]    1    2       

Answer 2

Here's a pretty basic Rcpp implementation:

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector MaxCol(Rcpp::NumericMatrix m) {
    R_xlen_t nr = m.nrow(), nc = m.ncol(), i = 0;
    Rcpp::NumericVector result(nr);

    for ( ; i < nr; i++) {
        double current = m(i, 0);
        R_xlen_t idx = 0, j = 1;
        for ( ; j < nc; j++) {
            if (m(i, j) > current) {
                current = m(i, j);
                idx = j;
            }
        }
        result[i] = idx + 1;
    }
    return result;
}

/*** R

microbenchmark::microbenchmark(
    "Rcpp" = MaxCol(Dummy), 
    "R" = max.col(Dummy, "first"),
    times = 200L
)
#Unit: milliseconds
# expr      min       lq     mean   median       uq      max neval
# Rcpp 221.7777 224.7442 242.0089 229.6407 239.6339 455.9549   200
# R    513.4391 524.7585 562.7465 539.4829 562.3732 944.7587   200

*/

---

I had to scale your sample data down by an order of magnitude since my laptop did not have enough memory, but the results should translate on your original sample data:

Dummy <- matrix(runif(50000000,0,3), ncol = 10000)
all.equal(MaxCol(Dummy), max.col(Dummy, "first"))
#[1] TRUE

---

This can be changed slightly to return the indices of the min and max in each row:

// [[Rcpp::export]]
Rcpp::NumericMatrix MinMaxCol(Rcpp::NumericMatrix m) {
    R_xlen_t nr = m.nrow(), nc = m.ncol(), i = 0;
    Rcpp::NumericMatrix result(nr, 2);

    for ( ; i < nr; i++) {
        double cmin = m(i, 0), cmax = m(i, 0);
        R_xlen_t min_idx = 0, max_idx = 0, j = 1;
        for ( ; j < nc; j++) {
            if (m(i, j) > cmax) {
                cmax = m(i, j);
                max_idx = j;
            }
            if (m(i, j) < cmin) {
                cmin = m(i, j);
                min_idx = j;
            }
        }
        result(i, 0) = min_idx + 1;
        result(i, 1) = max_idx + 1;
    }
    return result;
}

---