R - 如何从距离矩阵中获取匹配元素的行列下标

Question

我有一个整型向量 vec1 ，我正在使用 dist 函数生成一个距离矩阵。我想要获取特定数值元素在距离矩阵中的坐标（行和列）。实际上，我想要得到那些距离为d的一对元素。例如：

vec1 <- c(2,3,6,12,17)
distMatrix <- dist(vec1)

#   1  2  3  4
#2  1         
#3  4  3      
#4 10  9  6   
#5 15 14 11  5

假设我想找到向量中相距5个单位的一对元素。我需要获取距离矩阵的行坐标（coordinate1）和列坐标（coordinate2）。在这个示例中，我期望结果如下：

coord1  
# [1] 5
coord2
# [1] 4

我想知道是否有一种高效的方法来获取这些值，而不需要将dist对象转换为矩阵或循环遍历矩阵？

- DKangeyan

Answer 1

距离矩阵是以压缩格式存储的下三角矩阵，其中下三角矩阵按列存储为一维向量。 您可以通过以下方式进行检查：

str(distMatrix)
# Class 'dist'  atomic [1:10] 1 4 10 15 3 9 14 6 11 5
# ...

即使我们调用 dist(vec1, diag = TRUE, upper = TRUE)，向量仍然是相同的；只是打印风格发生了变化。也就是说，无论您如何调用 dist，您始终会得到一个向量。

本答案重点介绍如何在一维索引和二维索引之间进行转换，以便您可以在不使用 as.matrix 将其作为完整矩阵的情况下处理 "dist" 对象。如果您确实想将其转换为矩阵，请使用as.matrix on a distance object is extremely slow; how to make it faster?中定义的 dist2mat 函数。

---

---

R 函数

编写基于向量的 R 函数来进行这些索引转换很容易。我们只需要小心地处理 "out-of-bound" 索引，对于这种情况应该返回 NA。

## 2D index to 1D index
f <- function (i, j, dist_obj) {
  if (!inherits(dist_obj, "dist")) stop("please provide a 'dist' object")
  n <- attr(dist_obj, "Size")
  valid <- (i >= 1) & (j >= 1) & (i > j) & (i <= n) & (j <= n)
  k <- (2 * n - j) * (j - 1) / 2 + (i - j)
  k[!valid] <- NA_real_
  k
  }

## 1D index to 2D index
finv <- function (k, dist_obj) {
  if (!inherits(dist_obj, "dist")) stop("please provide a 'dist' object")
  n <- attr(dist_obj, "Size")
  valid <- (k >= 1) & (k <= n * (n - 1) / 2)
  k_valid <- k[valid]
  j <- rep.int(NA_real_, length(k))
  j[valid] <- floor(((2 * n + 1) - sqrt((2 * n - 1) ^ 2 - 8 * (k_valid - 1))) / 2)
  i <- j + k - (2 * n - j) * (j - 1) / 2
  cbind(i, j)
  }

这些函数非常节省内存使用，因为它们使用索引而不是矩阵。

---

将 finv 应用于你的问题

你可以使用

vec1 <- c(2,3,6,12,17)
distMatrix <- dist(vec1)

finv(which(distMatrix == 5), distMatrix)
#     i j
#[1,] 5 4

一般来说，距离矩阵包含浮点数。使用 == 判断两个浮点数是否相等是有风险的。请参考为什么这些数字不相等？获取更多信息和可能的策略。

---

使用 dist2mat 替代

使用给出的 dist2mat 函数替代 as.matrix，可以极大地提高速度。我们可以使用 which(, arr.ind = TRUE)。

library(Rcpp)
sourceCpp("dist2mat.cpp")
mat <- dist2mat(distMatrix, 128)
which(mat == 5, arr.ind = TRUE)
#  row col
#5   5   4
#4   4   5

附录：Markdown（需要MathJax支持）中的图片

---

## 2D index to 1D index

The lower triangular looks like this: $$\begin{pmatrix} 0 & 0 & \cdots & 0\\ \times & 0 & \cdots & 0\\ \times & \times & \cdots & 0\\ \vdots & \vdots & \ddots & 0\\ \times & \times & \cdots & 0\end{pmatrix}$$ If the matrix is $n \times n$, then there are $(n - 1)$ elements ("$\times$") in the 1st column, and $(n - j)$ elements in the j<sup>th</sup> column. Thus, for element $(i,\  j)$ (with $i > j$, $j < n$) in the lower triangular, there are $$(n - 1) + \cdots (n - (j - 1)) = \frac{(2n - j)(j - 1)}{2}$$ "$\times$" in the previous $(j - 1)$ columns, and it is the $(i - j)$<sup>th</sup> "$\times$" in the $j$<sup>th</sup> column. So it is the $$\left\{\frac{(2n - j)(j - 1)}{2} + (i - j)\right\}^{\textit{th}}$$ "$\times$" in the lower triangular.

----

## 1D index to 2D index

Now for the $k$<sup>th</sup> "$\times$" in the lower triangular, how can we find its matrix index $(i,\ j)$? We take two steps: 1> find $j$;  2> obtain $i$ from $k$ and $j$.

The first "$\times$" of the $j$<sup>th</sup> column, i.e., $(j + 1,\ j)$, is the $\left\{\frac{(2n - j)(j - 1)}{2} + 1\right\}^{\textit{th}}$ "$\times$" of the lower triangular, thus $j$ is the maximum value such that $\frac{(2n - j)(j - 1)}{2} + 1 \leq k$. This is equivalent to finding the max $j$ so that $$j^2 - (2n + 1)j + 2(k + n - 1) \geq 0.$$ The LHS is a quadratic polynomial, and it is easy to see that the solution is the integer no larger than its first root (i.e., the root on the left side): $$j = \left\lfloor\frac{(2n + 1) - \sqrt{(2n-1)^2 - 8(k-1)}}{2}\right\rfloor.$$ Then $i$ can be obtained from $$i = j + k - \left\{\frac{(2n - j)(j - 1)}{2}\right\}.$$

- Zheyuan Li

---

为了从“dist”矩阵中提取子对角线，可以使用函数f。然而，我特别编写了一个简化函数：在R中从距离矩阵中提取对角线。 - Zheyuan Li

Answer 2

如果向量不太大，最好的方法可能是将dist的输出包装到as.matrix中，并使用选项arr.ind=TRUE的which方法。在一个dist矩阵中检索索引号的标准方法的唯一缺点是增加了内存使用量，在传递给dist的非常大的向量的情况下，这可能变得重要。这是因为将dist返回的下三角矩阵转换为常规密集矩阵实际上会使存储的数据量翻倍。

另一种方法是将dist对象转换为列表，使dist的下三角矩阵中的每列表示列表的一个成员。然后可以将列表成员的索引编号和元素位置映射到密集N x N矩阵的列和行号，而不生成矩阵。

以下是这种基于列表的方法的一个可能的实现：

distToList <- function(x) {
  idx <- sum(seq(length(x) - 1)) - rev(cumsum(seq(length(x) - 1))) + 1
  listDist <- unname(split(dist(x), cumsum(seq_along(dist(x)) %in% idx)))
  # https://dev59.com/zWQo5IYBdhLWcg3wDbwJ#16358095
}
findDistPairs <- function(vec, theDist) {
  listDist <- distToList(vec)
  inList <- lapply(listDist, is.element, theDist)
  matchedCols <- which(sapply(inList, sum) > 0) 
  if (length(matchedCols) > 0) found <- TRUE else found <- FALSE
  if (found) {
    matchedRows <- sapply(matchedCols, function(x) which(inList[[x]]) + x ) 
  } else {matchedRows <- integer(length = 0)}
  matches <- cbind(col=rep(matchedCols, sapply(matchedRows,length)), 
                   row=unlist(matchedRows))
  return(matches)
}

vec1 <- c(2, 3, 6, 12, 17)
findDistPairs(vec1, 5)
#     col row
#[1,]   4   5

代码中可能不太清晰的部分涉及将列表中条目的位置映射到N x N矩阵的列/行值。虽然这些转换并不简单，但却很直观。

在代码中的注释中，我指出了一个在StackOverflow上的答案，该答案被用于将向量拆分为列表。循环（sapply、lapply）在性能方面应该不会有问题，因为它们的范围是O（N）级别的。此代码的内存使用情况主要由列表的存储确定。这种内存使用量应该与dist对象的相似，因为两个对象包含相同的数据。

dist对象在函数distToList（）中计算并转换为列表。由于dist计算在任何情况下都是必需的，所以在处理大型向量的情况下，此函数可能需要花费一定时间。如果目标是找到具有不同距离值的多个配对，则最好仅针对给定向量计算listDist一次，并将结果列表存储，例如存储在全局环境中。

长话短说，处理这类问题的常规方法简单快捷。

distMatrix <- as.matrix(dist(vec1)) * lower.tri(diag(vec1))
which(distMatrix == 5, arr.ind = TRUE)
#  row col
#5   5   4

我建议默认使用这种方法。在内存限制达到时，比如对于非常大的向量vec1的情况下，可能需要使用更复杂的解决方案。上面描述的基于列表的方法可以提供一种解决方法。

- RHertel

---

太棒了！我开始学习这个，发现它真的很有帮助。 - Geparada