Faster Hamming distance in R (2)

In my last post, I proposed a very fast and simple matrix multiplication approach for calculating the Hamming distance between all columns of a matrix X. This approach was restricted to binary matrices:

```
hamming_binary <- function(X, Y = NULL) {
if (is.null(Y)) {
D <- t(1 - X) %*% X
D + t(D)
} else {
t(1 - X) %*% Y + t(X) %*% (1 - Y)
}
}
```

Since it only relies on low-level matrix algebra routines, the above function is extremely fast, much faster for example than the hamming.distance function from the e1071 package, which is implemented as an R-level nested loop and is also restricted to binary matrices. However, what if we would like to calculate the Hamming distance between strings of characters, such as DNA sequences, which are composed of As, Cs, Gs and Ts? The nice thing is, that matrix multiplication can also be used for matrices with more than two possible values, and different types of elements. Furthermore, the approach below is generalized to calculating the Hamming distances between all pairs of columns of two matrices X and Y:

```
hamming <- function(X, Y = NULL) {
if (is.null(Y)) {
uniqs <- unique(as.vector(X))
U <- X == uniqs[1]
H <- t(U) %*% U
for ( uniq in uniqs[-1] ) {
U <- X == uniq
H <- H + t(U) %*% U
}
} else {
uniqs <- union(X, Y)
H <- t(X == uniqs[1]) %*% (Y == uniqs[1])
for ( uniq in uniqs[-1] ) {
H <- H + t(X == uniq) %*% (Y == uniq)
}
}
nrow(X) - H
}

```

The function above calculates the Hamming distance between all columns of a matrix X, or two matrices X and Y. Again matrix multiplication is used, this time for counting, between two columns x and y, the number of cases in which corresponding elements have the same value (e.g. A, C, G or T). This counting is done for each of the possible values individually, while iteratively adding the results. The end result of the iterative adding is the sum of all corresponding elements that are the same, i.e. the inverse of the Hamming distance. Therefore, the last step is to subtract this end result H from the maximum possible distance, which is the number of rows of matrix X.

Interestingly, such a "multi-level" Hamming distance can also be quite elegantly expressed as a sum of multiple binary Hamming distances:

```
hamming <- function(X, Y = NULL) {
if (is.null(Y)) {
uniqs <- unique(as.vector(X))
H <- hamming_binary(X == uniqs[1])
for ( uniq in uniqs[-1] ) {
H <- H + hamming_binary(X == uniq)
}
} else {
uniqs <- union(X, Y)
H <- hamming_binary(X == uniqs[1], Y == uniqs[1])
for ( uniq in uniqs[-1] ) {
H <- H + hamming_binary(X == uniq, Y == uniq)
}
}
H / 2
}

```

Note that at the end, we need to divide by two, because by summing the individual binary Hamming distances, we are counting each element twice.

Also note that the function is not dependent on the mode of the matrix. In other words, the matrix can consist of characters, integers, and any other type of element you can store in an R matrix.

15 thoughts on “Faster Hamming distance in R (2)”

1. Anonymous says:

It seems your “hamming” function only works on DNA sequences, but not on amino acid sequences where there are 20 different letters?

Like

1. Anonymous says:

Sorry, forgot to say that they were great posts!

Like

1. Thanks! Glad you like them.

For me, the hamming function above does in fact work for any type of element (and number of distinct values) you can store in a matrix. For example, try running the following code:

n <- 1000
m <- 100
X <- matrix(sample(LETTERS, n * m, replace = TRUE), nrow = n, ncol = m)
Y <- matrix(sample(LETTERS, n * m, replace = TRUE), nrow = n, ncol = m)
H <- hamming(X, Y)

# check a random example
i <- sample(m, 1)
j <- sample(m, 1)
sum(X[, i] != Y[, j]) == H[i, j]

Like

2. Anonymous says:

Very nice! My bad. Didn’t read the post carefully. My sequences in my input matrix are arranged by rows. After transposing it, everything works!

Thanks!

Like

3. chris2016 says:

Can you add an example to get weighted hamming? Where the weights are the number of levels per feature then each distance is divided by this?

Like

1. Assuming that by “feature” you mean “row” in the matrices X and Y above (not really clear from your question), feature-wise weights are implemented as follows:

hamming_binary <- function(X, w = rep(1, nrow(X))) {
w <- sqrt(w)
D <- t(w * (1 – X)) %*% (w * X)
D + t(D)
}

This is easily extended to the more complex hamming() function. I might write a separate post on this some time.

Like

1. chris2016 says:

Yeah that’s it. Thanks. Nice work.

Like

2. chris2016 says:

Sorry, but what if I want to get the feature weights for more levels, i.e. categories. I have from your other code this function I need to adapt for my package, so how to change this?

hamming <- function(X) {
uniqs <- unique(as.vector(X))
U <- X == uniqs[1]
H <- t(U) %*% U
for ( uniq in uniqs[-1] ) {
U <- X == uniq
H <- H + t(U) %*% U
}
nrow(X) – H
}

Like

1. chris2016 says:

Ok this is good. I’ll let you know when this progresses and certainly can cite blog.

Like