R - 如何从距离矩阵中获取匹配元素的行和列下标
R - How to get row & column subscripts of matched elements from a distance matrix
我有一个整数向量 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 对象转换为矩阵或循环遍历矩阵的有效方法来获取这些值？
距离矩阵是压缩格式的下三角矩阵，其中下三角按列存储为一维向量。您可以通过
查看
str(distMatrix)
# Class 'dist' atomic [1:10] 1 4 10 15 3 9 14 6 11 5
# ...
即使我们调用dist(vec1, diag = TRUE, upper = TRUE)，向量还是一样的；只有打印样式发生变化。即无论怎么调用dist，得到的总是一个vector
这个答案着重于如何在 1D 和 2D 索引之间进行转换，这样您就可以使用 "dist" 对象，而无需先使用 as.matrix 使其成为完整矩阵。如果你确实想让它成为一个矩阵，使用 .
中定义的 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
一般来说，距离矩阵包含浮点数。用==判断两个浮点数是否相等是有风险的。阅读 Why are these numbers not equal? 了解更多可能的策略。
替代dist2mat
利用中给出的dist2mat函数，我们可以使用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\}.$$
如果向量不是太大，最好的方法可能是将 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)))
#
}
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 矩阵的列/行值的映射。这些转换虽然不简单，但很简单。
在代码中的评论中，我指出了 Whosebug 上的一个答案，该答案已在此处用于将矢量拆分为列表。循环 (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 的情况下，可能需要更复杂的解决方案。然后，上述基于列表的方法可以提供补救措施。
我有一个整数向量 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 对象转换为矩阵或循环遍历矩阵的有效方法来获取这些值？
距离矩阵是压缩格式的下三角矩阵，其中下三角按列存储为一维向量。您可以通过
查看str(distMatrix)
# Class 'dist' atomic [1:10] 1 4 10 15 3 9 14 6 11 5
# ...
即使我们调用dist(vec1, diag = TRUE, upper = TRUE)，向量还是一样的；只有打印样式发生变化。即无论怎么调用dist，得到的总是一个vector
这个答案着重于如何在 1D 和 2D 索引之间进行转换，这样您就可以使用 "dist" 对象，而无需先使用 as.matrix 使其成为完整矩阵。如果你确实想让它成为一个矩阵，使用
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
一般来说，距离矩阵包含浮点数。用==判断两个浮点数是否相等是有风险的。阅读 Why are these numbers not equal? 了解更多可能的策略。
替代dist2mat
利用dist2mat函数，我们可以使用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\}.$$
如果向量不是太大，最好的方法可能是将 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)))
#
}
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 矩阵的列/行值的映射。这些转换虽然不简单，但很简单。
在代码中的评论中，我指出了 Whosebug 上的一个答案，该答案已在此处用于将矢量拆分为列表。循环 (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 的情况下，可能需要更复杂的解决方案。然后，上述基于列表的方法可以提供补救措施。