R - 如何从距离矩阵中得到匹配元素的行和列下标

Question

我有一个整数向量vec1，我使用dist函数生成一个远距离矩阵。我想获得距离矩阵中某个值元素的坐标（行和列）。本质上，我希望得到一对相距甚远的元素。例如：

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个单位的向量中的元素对。我想获得坐标1，它们是行和坐标2，它们是距离矩阵的列。在这个玩具例子，我希望

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

我想知道是否有一个有效的方式来获取这些值不涉及dist对象转换为一个矩阵或循环通过矩阵？

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，您总是会获得一维数组。

假设一个完整的下三角是n-by-n，那么它的第(i,j)个元素将被映射到打包的1D数组中的第(j - 1) * (2 * n - 2 - j)/2 + (i - 1)个元素。我们可以定义一个指数变换函数：

## `i` and `j` can both be vector input, but they must have the same length 
f <- function (i, j, n) { 
    ifelse((i > j) & (j <= n), (j - 1) * (2 * n - 2 - j)/2 + (i - 1), NA_real_) 
    } 

在另一方面，如果我们知道包装的数组中元素的位置，说k，我们可以通过一个稍微复杂的功能找到(i,j)：

## `k` can be a vector input 
finv <- function (k, n) { 
    ## starting position for each column 
    ptr_all_cols <- f(2:n, 1:(n - 1), n) 
    ## maximum valid `k` 
    k_max <- n * (n - 1)/2 
    ## `finv` operation on a scalar `k` 
    scaler_finv <- function (k) { 
    if (k > k_max) return(c(i = NA_real_, j = NA_real_)) 
    j <- sum(ptr_all_cols <= k) ## get column index j 
    i <- k - ptr_all_cols[j] + j + 1 ## get row index i 
    c(i = i, j = j) 
    } 
    ## "vectorization" 
    do.call(rbind, lapply(k, scaler_finv)) 
    } 

这些转换函数在内存使用上非常便宜，因为它们使用索引而不是矩阵。

---

基于变换函数finv

随着finv有效的解决方案，它是晚饭有效地找到所需的元素。对于你的玩具例如，你可以使用

## the first `5` is the value to be matched; the second is matrix dimension 
finv(which(distMatrix == 5), 5) 
#  i j 
#[1,] 5 4 

---

注意

一般来说，距离矩阵包含浮点数。使用==来判断两个浮点数是否相等是相当危险的。阅读Why are these numbers not equal?了解更多和可能的策略。

---

替代

有由@RHertel提出一个方便的答案。那些拥有10,000声誉仍然能够看到它：

mat <- stats:::as.matrix.dist(dist(vec1)) * lower.tri(diag(vec1)) 
which(mat == 5, arr.ind = TRUE) 

另一种方式把第一行是

mat <- matrix(0, n, n); mat[lower.tri(mat)] <- distMatrix 

无论哪种方式，将花费更多的内存矩阵过程中存储了许多n-by-n矩阵操作（虽然后者相对便宜）。当vec1很长时，内存问题可能是一个瓶颈。

---

其它

f和finv可能是广义上非常有用的功能，至少它可以帮助理解全格式和压缩格式之间的指标怎么可以相互转化。

以下两个函数仅用于演示目的，它还检查f和finv的正确性。

## a function to verbose `f` transform, primarily used to check the correctness of `f` 
verbose_f <- function (n) { 
    i <- rep(seq_len(n), times = n) 
    j <- rep(seq_len(n), each = n) 
    matrix(f(i, j, n), n) 
    } 

## a function to verbose `finv` transform, primarily used to check the correctness of `finv` 
verbose_finv <- function (k, n) cbind(k = k, finv(k, n)) 

我们以n = 5为例。

verbose_f(5) 

#  [,1] [,2] [,3] [,4] [,5] 
#[1,] NA NA NA NA NA 
#[2,] 1 NA NA NA NA 
#[3,] 2 5 NA NA NA 
#[4,] 3 6 8 NA NA 
#[5,] 4 7 9 10 NA 

verbose_finv(1:15,5) 

#  k i j 
# [1,] 1 2 1 
# [2,] 2 3 1 
# [3,] 3 4 1 
# [4,] 4 5 1 
# [5,] 5 3 2 
# [6,] 6 4 2 
# [7,] 7 5 2 
# [8,] 8 4 3 
# [9,] 9 5 3 
#[10,] 10 5 4 
#[11,] 11 NA NA 
#[12,] 12 NA NA 
#[13,] 13 NA NA 
#[14,] 14 NA NA 
#[15,] 15 NA NA 

在这两种情况下，NA暗示 “下标越界”。

Answer 2

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

另一种方法是将dist对象转换为列表，使得dist的下三角矩阵中的每列代表列表的一个成员。然后可以将列表成员的索引号和列表成员中的元素的位置映射到密集的N×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))) 
    # http://stackoverflow.com/a/16358095/4770166 
} 
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×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的情况下，可能需要更复杂的解决方案。然后上述的基于列表的方法可以提供补救措施。