计算两个整数 matrices/data 帧的所有行之间的成对汉明距离

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 的列，给出第二个最小距离；
距离第二小的列。

谢谢。

Answer 1

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

正如我在评论中所说，我们可以这样做：

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

计算两个等长整数向量x和y之间的汉明距离。这仅使用 R 基础，但比 e1071::hamming.distance、 更有效，因为它是矢量化的！

对于您 post 中的示例 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)))
    }
  }

注意：

hmd 执行计算按列。它被设计为 CPU 缓存友好 。这样，如果我们想做一些按行的计算，我们应该先转置矩阵；
这里没有明显的循环；相反，我们使用 tapply().

两个matrices/data帧之间的快速汉明距离计算

这就是你想要的。下面的函数 foo 接受两个数据框或矩阵 df1 和 df2，计算 df1 和 df2 的每一行之间的距离。 argument p 是一个整数，表示要保留多少个结果。 p = 3 将在 df1.

中保持最小的 3 个行 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)
  }

注意：

开头做换位，根据前面的原因；
这里使用了一个for循环。但这实际上是有效的，因为在每次迭代中都进行了大量计算。它也比使用 *apply 系列更优雅，因为我们要求多个输出（行 ID id 和距离 d）。

实验

这部分使用小型数据集来 test/demonstrate 我们的函数。

一些玩具数据：

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，您知道该怎么做，对吗？

Answer 2

请不要惊讶我为什么要选另一个部分。这部分给出了一些相关的东西。 这不是OP所要求的，但可能会对任何读者有所帮助。

一般汉明距离计算

在上一个答案中，我从一个函数 hmd0 开始，该函数计算两个相同长度的整数向量之间的汉明距离。这意味着如果我们有 2 个整数向量：

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。在 my answer to this question 中，我演示了使用编译代码的想法。这当然需要写一个C版的hamming.distance，这里就不展示了

Answer 3

这是一个仅使用基础 R 的替代解决方案，应该非常快，尤其是当您的 df1 和 df2 有很多行时。这样做的主要原因是它没有使用 any R 级循环来计算汉明距离，例如 for 循环、while 循环或 *apply 函数。相反，它使用 matrix multiplication for computing the Hamming distance. In R, this is much faster than any approach using R-level looping. Also note that using an *apply function will not necessarily make your code any faster than using a for loop. Two other efficiency-related features of this approach are: (1) It uses partial sorting 为 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 秒。

计算两个整数 matrices/data 帧的所有行之间的成对汉明距离

Computing pairwise Hamming distance between all rows of two integer matrices/data frames

r

apply

hamming-distance

sapply

tapply