# Computing pairwise distance covariances for a symmetric matrix

I wrote a non-vectorized, non-parallel wrapper function for computing the pairwise distance covariances (Wikipedia) (original paper) in a data frame:

library(energy)

pairwise_dcor <- function (DF, progress = FALSE, parallel = FALSE) {
stopifnot(is.data.frame(DF))
v <- names(DF)
n <- length(v)
corr <- matrix(1, n, n, dimnames = list(v, v))
for (i in seq_len(n)) {
for (j in seq_len(i - 1)) {
cc <- complete.cases(DF[[i]], DF[[j]])
xi <- DF[cc, i]
xj <- DF[cc, j]
if (is.factor(xi) || is.character(xi)) {
xi <- dist(model.matrix(~ xi - 1), "binary")
}
if (is.factor(xj) || is.character(xj)) {
xj <- dist(model.matrix(~ xj - 1), "binary")
}
corr[i, j] <- dcor(xi, xj)
corr[j, i] <- corr[i, j]
}
}
corr
}


It's very slow:

# generate data that's similar to mine
test_data <- do.call(rbind, replicate(20, iris, simplify = FALSE))
test_data <- do.call(cbind, replicate(5, test_data, simplify = FALSE))

# run it (but i don't recommend it)
st <- system.time(
dc <- pairwise_dcor(test_data)
)
st
#     user   system  elapsed
#  955.179  102.529 1075.781


This test data has a run time of several minutes. I need to repeat this operation several (10-20) times, and a few times on a much larger (like 3,000 x 100) data set, with missing values in most or all of the columns. I'd rather not spend most of the afternoon waiting for this code to run, since I'd like to be able to run it on-the-fly as I update my analysis.

It seems like it would be easy to improve its performance by parallelizing it. I know about the doParallel and foreach packages, but I'm not sure what the best way to use it would be and I've had mixed results with it in the past, with more cumbersome, less portable code and dubious performance improvements. For what it's worth, I'd like to run this on a quad-core MacBook Pro with 16 GB of RAM.

Also, am I missing out on other optimizations? For instance, the dist(model.matrix(~ xi - 1), "binary") line seems absurdly roundabout to me, but my brute-force reimplementation was substantially slower:

binary_jaccard <- function(x) {
n <- length(x)
d <- matrix(0, n, n)
for (i in seq_len(n)) {
for (j in seq_len(i - 1)) {
d[i, j] <- as.numeric(x[i] != x[j])
d[j, i] <- d[i, j]
}
}
as.dist(d)
}

silly_version <- as.matrix(dist(model.matrix(~iris$Species - 1), "binary")) slow_version <- binary_jaccard(iris$Species, FALSE)
all(silly_version == slow_version)
# TRUE

system.time(as.matrix(dist(model.matrix(~iris$Species - 1), "binary"))) # user system elapsed # 0.002 0.000 0.003 system.time(binary_jaccard(iris$Species, FALSE))
#    user  system elapsed
#   0.988   0.071   1.061


Could I see improvements by using Rcpp?

Alternatively, would I be able to make significant improvements by switching to Python+NumPy+Pandas?

• This should give you ideas: stackoverflow.com/a/27943600/1201032, it only uses fast R functions written in C. Also, the converting of your data into numeric is something you only need to do once, at the very beginning of your function. – flodel Mar 11 '15 at 10:33
• @flodel just to be clear, "distance correlation" is not Pearson correlation. Are you suggesting that I eschew the dcor function entirely (written in C btw)? – shadowtalker Mar 11 '15 at 13:40
• Ok... Reminded me of something else then: stackoverflow.com/questions/13973185/…. Not a full working solution yet, just an idea. – flodel Mar 12 '15 at 0:12
• @flodel that's a lot more helpful, thanks. Unfortunately my data is riddled with missing values, so I need to recompute A and B every iteration in order to mimic the "pairwise.complete" option in cor. – shadowtalker Mar 12 '15 at 2:45
• @flodel I took that advice and made a follow-up question that is also somewhat more coherent – shadowtalker Mar 12 '15 at 7:07