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的情況下，可能需要更復雜的解決方案。然後上述的基於列表的方法可以提供補救措施。