计算两个整数矩阵/数据框中所有行之间的汉明距离

Question

我有两个数据框，其中df1是参考数据，df2是新数据。对于df2中的每一行，我需要找到与df1中最匹配（和次匹配）的行，以汉明距离为依据。

我使用e1071包来计算汉明距离。例如，可以通过以下方式计算两个向量x和y之间的汉明距离：

x <- c(356739, 324074, 904133, 1025460, 433677, 110525, 576942, 526518, 299386,
       92497, 977385, 27563, 429551, 307757, 267970, 181157, 3796, 679012, 711274,
       24197, 610187, 402471, 157122, 866381, 582868, 878)

y <- c(356739, 324042, 904133, 959893, 433677, 110269, 576942, 2230, 267130,
       92496, 960747, 28587, 429551, 438825, 267970, 181157, 36564, 677220,
       711274, 24485, 610187, 404519, 157122, 866413, 718036, 876)

xm <- sapply(x, intToBits)
ym <- sapply(y, intToBits)

distance <- sum(sapply(1:ncol(xm), function(i) hamming.distance(xm[,i], ym[,i])))

结果距离为25。但我需要对df1和df2的所有行都执行此操作。一个简单的方法需要双重循环嵌套，看起来非常慢。

有什么更有效的方法吗？最终我需要附加到df2：

• 一列是df1中给出最小距离的行id；
• 一列是最小距离；
• 一列是df1中给出第二小距离的行id；
• 一列是第二小距离。

谢谢。

- alaj

---

应该能够使用 apply 和 match 来完成它。 - Hack-R

Answer 1

快速计算等长整数向量之间的汉明距离

正如我在评论中所说，我们可以进行以下操作：

hmd0 <- function(x,y) sum(as.logical(xor(intToBits(x),intToBits(y))))

计算两个长度相等的整数向量之间的海明距离。这仅使用R基础知识，但比“e1071 :: hamming.distance”更有效，因为它是矢量化的！对于您发布的示例“x”和“y”，这将给出25。（如果我们想要成对的海明距离，我的另一个答案将展示我们应该怎么做。）

快速计算矩阵和向量之间的海明距离

如果我们想要计算单个“y”和多个“x”之间的海明距离，即向量和矩阵之间的海明距离，我们可以使用以下函数。

hmd <- function(x,y) {
  rawx <- intToBits(x)
  rawy <- intToBits(y)
  nx <- length(rawx)
  ny <- length(rawy)
  if (nx == ny) {
    ## quick return
    return (sum(as.logical(xor(rawx,rawy))))
    } else if (nx < ny) {
    ## pivoting
    tmp <- rawx; rawx <- rawy; rawy <- tmp
    tmp <- nx; nx <- ny; ny <- tmp
    }
  if (nx %% ny) stop("unconformable length!") else {
    nc <- nx / ny  ## number of cycles
    return(unname(tapply(as.logical(xor(rawx,rawy)), rep(1:nc, each=ny), sum)))
    }
  }

请注意:

1. hmd是按列进行计算的。它被设计成“CPU缓存友好”的方式。因此，如果我们想要进行一些按行计算，我们应该先转置矩阵; 2. 这里没有明显的循环;而是使用了tapply()。

快速计算两个矩阵/数据帧之间的汉明距离

这就是您想要的。以下函数foo接受两个数据帧或矩阵df1和df2，计算df1和df2每一行之间的距离。参数p是一个整数，显示您想要保留多少结果。p = 3将保留最小的3个距离及其在df1中的行ID。

foo <- function(df1, df2, p) {
  ## check p
  if (p > nrow(df2)) p <- nrow(df2)
  ## transpose for CPU cache friendly code
  xt <- t(as.matrix(df1))
  yt <- t(as.matrix(df2))
  ## after transpose, we compute hamming distance column by column
  ## a for loop is decent; no performance gain from apply family
  n <- ncol(yt)
  id <- integer(n * p)
  d <- numeric(n * p)
  k <- 1:p
  for (i in 1:n) {
    distance <- hmd(xt, yt[,i])
    minp <- order(distance)[1:p]
    id[k] <- minp
    d[k] <- distance[minp]
    k <- k + p
    }
  ## recode "id" and "d" into data frame and return
  id <- as.data.frame(matrix(id, ncol = p, byrow = TRUE))
  colnames(id) <- paste0("min.", 1:p)
  d <- as.data.frame(matrix(d, ncol = p, byrow = TRUE))
  colnames(d) <- paste0("mindist.", 1:p)
  list(id = id, d = d)
  }

请注意:

1. 在开始时进行置换，根据之前的原因； 2. 在这里使用了一个for循环。但是，这实际上是有效的，因为每个迭代中都进行了相当多的计算。而且，与使用*apply家族相比，这更优雅，因为我们要求多个输出(行id id和距离d)。

---

实验

本部分使用小数据集来测试/演示我们的函数。

一些玩具数据:

set.seed(0)
df1 <- as.data.frame(matrix(sample(1:10), ncol = 2))  ## 5 rows 2 cols
df2 <- as.data.frame(matrix(sample(1:6), ncol = 2))  ## 3 rows 2 cols

首先测试hmd(需要转位)：

hmd(t(as.matrix(df1)), df2[1, ])  ## df1 & first row of df2
# [1] 2 4 6 2 4

测试foo：

foo(df1, df2, p = 2)

# $id
#   min1 min2
# 1    1    4
# 2    2    3
# 3    5    2

# $d
#   mindist.1 mindist.2
# 1         2         2
# 2         1         3
# 3         1         3

如果你想在df2中添加一些列，你知道该怎么做，对吗？

- Zheyuan Li

---

非常感谢。你所做的非常清晰明了。我发现foo函数存在一个问题，就是在代码末尾硬编码了ncol为3。我认为你应该将其设置为p。 - alaj

当然可以。再次感谢。我还在尝试弄清楚如何集成另外两个数字：df2和最低距离df1中设置为1的位数。我需要一个新函数来完成这个任务吗？还是可以将其集成到hmd函数中？您有什么指导可以提供吗？ - alaj

谢谢。我已经创建了一个名为“计算两个数据帧之间汉明距离匹配行中被设置为1的位数”的新帖子。 - alaj

Answer 2

请不要惊讶为什么我另起一段。这部分提供了相关信息。虽然不是 OP 所要求的，但可能对任何读者有所帮助。

---

通用海明距离计算

在之前的回答中，我从一个函数 hmd0 开始，它可以计算两个相同长度的整数向量之间的海明距离。这意味着如果我们有两个整数向量：

set.seed(0)
x <- sample(1:100, 6)
y <- sample(1:100, 6)

我们最终会得到一个标量：

hmd0(x,y)
# 13

如果我们想要计算两个向量的逐对汉明距离，该怎么办呢？

实际上，我们的函数hmd只需要进行简单的修改就可以实现：

hamming.distance <- function(x, y, pairwise = TRUE) {
  nx <- length(x)
  ny <- length(y)
  rawx <- intToBits(x)
  rawy <- intToBits(y)
  if (nx == 1 && ny == 1) return(sum(as.logical(xor(intToBits(x),intToBits(y)))))
  if (nx < ny) {
    ## pivoting
    tmp <- rawx; rawx <- rawy; rawy <- tmp
    tmp <- nx; nx <- ny; ny <- tmp
    }
  if (nx %% ny) stop("unconformable length!") else {
    bits <- length(intToBits(0)) ## 32-bit or 64 bit?
    result <- unname(tapply(as.logical(xor(rawx,rawy)), rep(1:ny, each = bits), sum))
    }
  if (pairwise) result else sum(result)
  }

现在

hamming.distance(x, y, pairwise = TRUE)
# [1] 0 3 3 2 5 0
hamming.distance(x, y, pairwise = FALSE)
# [1] 13

---

海明距离矩阵

如果我们想要计算海明距离矩阵，例如：

set.seed(1)
x <- sample(1:100, 5)
y <- sample(1:100, 7)

x和y之间的距离矩阵为:

outer(x, y, hamming.distance)  ## pairwise argument has no effect here

#      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
# [1,]    2    3    4    3    4    4    2
# [2,]    7    6    3    4    3    3    3
# [3,]    4    5    4    3    6    4    2
# [4,]    2    3    2    5    6    4    2
# [5,]    4    3    4    3    2    0    2

我们还可以这样做：

outer(x, x, hamming.distance)

#     [,1] [,2] [,3] [,4] [,5]
# [1,]    0    5    2    2    4
# [2,]    5    0    3    5    3
# [3,]    2    3    0    2    4
# [4,]    2    5    2    0    4
# [5,]    4    3    4    4    0

在后一种情况下，我们最终得到一个对角线上为0的对称矩阵。在这里使用outer是低效的，但它仍然比编写R循环更有效率。由于我们的hamming.distance是用R代码编写的，所以我会继续使用outer。在我的回答和这个问题中，我演示了使用已编译代码的思路。当然，这需要编写hamming.distance的C版本，但我不会在这里展示它。

- Zheyuan Li

Answer 3

这是一种使用纯基础 R 的替代方案，特别适用于 df1 和 df2 行数较多的情况下，且速度非常快。主要原因是它不使用任何 R 级别的循环来计算汉明距离，例如 for 循环、while 循环或 *apply 函数。相反，它使用 矩阵乘法计算汉明距离。在 R 中，这比使用 R 级别循环的任何方法都要快得多。还要注意的是，使用 *apply 函数并不一定比使用 for 循环更快。此方法的另外两个与效率相关的特点是：(1) 它使用 部分排序 来查找 df2 中每行的最佳两个匹配项，以及 (2) 它将 df1 的整个位表示存储在一个矩阵中（df2 也是如此），并且可以在一个步骤中完成，而不使用任何 R 级别的循环。

执行所有工作的函数：

# INPUT:       
# X corresponds to your entire df1, but is a matrix
# Y corresponds to your entire df2, but is a matrix
# OUTPUT:
# Matrix with four columns corresponding to the values 
# that you specified in your question
fun <- function(X, Y) {

  # Convert integers to bits 
  X <- intToBits(t(X))
  # Reshape into matrix
  dim(X) <- c(ncols * 32, nrows)

  # Convert integers to bits
  Y <- intToBits(t(Y))
  # Reshape into matrix
  dim(Y) <- c(ncols * 32, nrows)

  # Calculate pairwise hamming distances using matrix 
  # multiplication. 
  # Columns of H index into Y; rows index into X.
  # The code for the hamming() function was retrieved
  # from this page:
  # https://johanndejong.wordpress.com/2015/10/02/faster-hamming-distance-in-r-2/
  H <- hamming(X, Y)

  # Now, for each row in Y, find the two best matches 
  # in X. In other words: for each column in H, find 
  # the two smallest values and their row indices.
  t(apply(H, 2, function(h) {
    mindists <- sort(h, partial = 1:2)
    c(
      ind1 = which(h == mindists[1])[1],
      val1 = mindists[1],
      hmd2 = which(h == mindists[2])[1],
      val2 = mindists[2]
    )
  }))
}

调用函数处理一些随机数据：

# Generate some random test data with no. of columns 
# corresponding to your data
nrows <- 1000
ncols <- 26 

# X corresponds to your df1
X <- matrix(
  sample(1e6, nrows * ncols, replace = TRUE), 
  nrow = nrows, 
  ncol = ncols
)

# Y corresponds to your df2
Y <- matrix(
  sample(1e6, nrows * ncols, replace = TRUE), 
  nrow = nrows, 
  ncol = ncols
)

res <- fun(X, Y)

上述示例中，X（df1）和Y（df2）均包含1000行，在我的笔记本电脑上运行大约需要1.1-1.2秒。

- Johann de Jong