3

I am looking for a more efficient way to get the distance matrix in terms of Hamming distance.

Backgrounds

I know there is a function hamming.distance() from package e1071 to compute the distance matrix, but I suspect it might be very slow when involving a large matrix with many rows, since it applied nested for loops for computation.

So far I have a faster way (see methodB) in the code below. However, it is only suitable for in the binary domain, i.e., {0,1}^n. However, it is unavailable when encountering domains consisting of more than 2 elements, i.e., {0,1,2,...,K-1}^n. In this sense, methodB is not for generic hamming distance.

Objective

My objective is to find a approach having the following features:

  • composed by functions only from base R (not using Rcpp to rewrite function for speeding up)
  • faster than my approach methodB() for the special case k=2
  • can be generalized for any positive integer k
  • outperform the speed of hamming.distance() from package e1071

My code

library(e1071)
# vector length, i.e., number of matrix
n <- 7
# number of elements to consist of domain {0,1,...,k-1}^n
k <- 2
# matrix for computing hamming distances by rows
m <- as.matrix(do.call(expand.grid,replicate(n,list(0:k-1))))

# applying `hamming.distance()` from package "e1071", which is generic so it is available for any positive integer `k`
methodA <- function(M) hamming.distance(M)
# my customized method from base R function `dist()`, which is not available for cases `k >= 2`
methodB <- function(M) as.matrix(round(dist(M,upper = T,diag = T)**2))

and the benchmark gives

microbenchmark::microbenchmark(
  methodA(m),
  methodB(m),
  unit = "relative",
  check = "equivalent",
  times = 50
)

Unit: relative
       expr      min       lq   mean   median       uq      max neval
 methodA(m) 33.45844 33.81716 33.963 34.30313 34.92493 14.92111    50
 methodB(m)  1.00000  1.00000  1.000  1.00000  1.00000  1.00000    50

Appreciated in advance!

ThomasIsCoding
  • 96,636
  • 9
  • 24
  • 81

3 Answers3

3

I found this blog which has four posts on calculcating Hamming Matrixes. I don't want to claim any fame for it, but maybe have a look at it. https://johanndejong.wordpress.com/2015/10/02/faster-hamming-distance-in-r-2/

hamming <- function(X) {
  D <- (1 - X) %*% t(X)
  D + t(D)
}

> microbenchmark::microbenchmark(
+   methodB(m),
+   hamming(m),
+   unit = "relative",
+   times = 50
+ )
Unit: relative
       expr    min       lq     mean   median       uq      max neval
 methodB(m) 1.0000 1.000000 1.000000 1.000000 1.000000 1.000000    50
 hamming(m) 1.2502 1.299844 1.436486 1.301461 1.302033 4.607748    50

PS: I don't have enough reputation to just leave this as a comment.

arohland
  • 96
  • 6
2
methodM <- function(x) {
  xt <- t(x)
  sapply(1:nrow(x), function(y) colSums(xt != xt[, y]))
}
microbenchmark::microbenchmark(
  methodB(m), methodM(m),
  unit = "relative", check = "equivalent", times = 50
)
# Unit: relative
#       expr  min       lq     mean   median       uq      max neval cld
# methodB(m) 1.00 1.000000 1.000000 1.000000 1.000000 1.000000    50  a 
# methodM(m) 1.25 1.224827 1.359573 1.219507 1.292463 4.550159    50   b
minem
  • 3,640
  • 2
  • 15
  • 29
  • Thanks for the answer! Yes, your method is fast and can be extended to cases where `k>2`. Do you have any idea that is faster than my `methodB()` for the special case `k = 2`? – ThomasIsCoding Jan 09 '20 at 13:34
  • @ThomasIsCoding it isn't faster probably because doesn't use low level, optimized function `dist`. Also this calculates all distances, ideal we would want to use better algorithm that calculates the values only for upper(lower) triangular part. – minem Jan 14 '20 at 08:13
  • yes, I see, maybe what we have is already the speed limit of using base R – ThomasIsCoding Jan 14 '20 at 08:18
-1

Did you try using Rcpp? I had a very similar problem! Please see the answer here: https://stackoverflow.com/a/60067825/3237589

Sudaraka
  • 125
  • 7