4
\$\begingroup\$

I would like to compute the number of concordant and discordant pairs of two large vectors. This function is a nice try, and works quite efficiently.

pairs <- function(x,y){
    n <- length(x)
    ix <- order(x)
    x <- x[ix]
    y <- y[ix]
    Nc <- sum(sapply(1:(n-1),function(i) sum(x[i]<x[(i+1):n] & y[i]<y[(i+1):n])))
    Nd <- sum(sapply(1:(n-1),function(i) sum(x[i]<x[(i+1):n] & y[i]>y[(i+1):n])))
    return(list(Nc=Nc,Nd=Nd))
}

x <- runif(10000)
y <- runif(10000)
system.time(pairs(x,y))

Do you have any idea on how it would be possible to boost the function?

\$\endgroup\$
4
  • \$\begingroup\$ Check ConDis.matrix function from asbio package \$\endgroup\$
    – m-dz
    Jul 6 '17 at 15:49
  • 2
    \$\begingroup\$ You've already sorted x right? So you don't need to calculate x[i]<x[(i+1):n]? Removing this calculation will half your computation time (through my testing). \$\endgroup\$
    – Chi Pak
    Jul 6 '17 at 15:50
  • 1
    \$\begingroup\$ It's a clever idea, and will work with the aforementioned example, but let's imagine we have some ties in the data, x <- sample(1:10, 1000, replace = T);y <- sample(1:10, 1000, replace = T). In this case the number of concordant and discordant pairs would be not correctly counted. \$\endgroup\$
    – And
    Jul 6 '17 at 16:11
  • \$\begingroup\$ @And Removing dupes for that matter? \$\endgroup\$
    – M--
    Jul 7 '17 at 19:58
2
\$\begingroup\$

These sorts of operations are tough to vectorize because you need to compare every element to all later elements. If the counting operation is a performance bottleneck in your code, then one possibility would be to implement it in C++ using the Rcpp package:

library(Rcpp)
cppFunction(
  "IntegerVector CDcount(NumericVector x, NumericVector y) {
    IntegerVector counts(2, 0);
    int n = x.size();
    for (int i=0; i < n; ++i) {
      for (int j=i+1; j < n; ++j) {
        counts[0] += (x[i] < x[j]) && (y[i] < y[j]);
        counts[1] += (x[i] < x[j]) && (y[i] > y[j]);
      }
    }
    return counts;
  }")

pairs2 <- function(x, y) {
  n <- length(x)
  ix <- order(x)
  x <- x[ix]
  y <- y[ix]
  counts <- CDcount(x, y)
  return(list(Nc=counts[1], Nd=counts[2]))
}

For your test data, this results in a 50 times speedup on my system.

x <- runif(10000)
y <- runif(10000)
system.time(pairs(x,y))
#    user  system elapsed 
#   4.892   1.447   6.584 
system.time(pairs2(x, y))
#    user  system elapsed 
#   0.119   0.001   0.120 
identical(pairs(x, y), pairs2(x, y))
# [1] TRUE

The Rcpp package takes some time to compile the CDcount function, so this option is probably only helpful if you are operating on very large vectors or if you are repeating the operation many times.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.