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

下三角が列によって1次元ベクトルとして格納され、下三角行列インデックス変換

距離マトリックスがパックフォーマットの下三角行列であり、のパックストレージ。あなたは私たちがdist(vec1, diag = TRUE, upper = TRUE)を呼び出す場合でも、結果は印刷スタイルの変更点を除いて、まだ同じであることを

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

注意を経由して、これを確認することができます。要約すると、distとどう呼んでも、常に1D配列が得られます。

はその(i,j)番目の要素が充填された1次元アレイで(j - 1) * (2 * n - 2 - j)/2 + (i - 1)番目の要素にマッピングされ、完全下三角はn-by-nであると仮定する。私たちはパック配列内の要素の位置を知っていれば、我々はもう少し複雑な関数で(i,j)を見つけることができ、kを言って、一方

## `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` 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 

---

注意

を使用することができ、距離行列は、浮動小数点数が含まれています。 2つの浮動小数点数が等しいかどうかを判断するのに、==を使うのはかなり危険です。より多くの可能な戦略については、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) 

もう一つの方法は、いずれかの方法では、マトリックス中にn-by-n行列の数を格納するため、より多くのメモリの費用がかかります

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

です（後者は比較的安価ですが）。 vec1が長いとメモリの問題がボトルネックになる可能性があります。

---

その他

機能fとfinvは、少なくとも、それは完全なフォーマットとパック形式の間でインデックスが相互変換することができますどのように理解するのに役立ち、広い意味ではかなり役に立つかもしれません。

次の2つの機能は、デモ目的のためのもので、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は「subscript out of bound」を意味します。

Answer 2

ベクトルが大きすぎない場合は、distの出力をas.matrixにラップし、arr.ind=TRUEの場合はwhichを使用することをお勧めします。 dist行列内のインデックス番号を検索するこの標準的な方法の唯一の欠点は、メモリ使用量の増加です。これは、非常に大きなベクトルがdistに渡された場合に重要になる場合があります。これは、distによって返された下三角行列の通常の密な行列への変換が、格納されたデータの量を事実上2倍にするためです。

代替案は、下位三角行列distの各列がリストの1メンバーを表すように、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の回答を指摘しました。ループ（サプリー、ラップリー）は、その範囲がオーダーO（N）であるため性能面で問題がないはずです。このコードのメモリ使用量は、主にリストの格納場所によって決まります。両方のオブジェクトに同じデータが含まれているため、この量のメモリはdistオブジェクトのメモリと似ていなければなりません。

distオブジェクトが計算され、関数distToList()のリストに変換されます。いずれにしても必要なdist計算のために、この関数は大きなベクトルの場合には時間がかかることがあります。異なる距離値を有するいくつかのペアを見つけることを目標とする場合、与えられたベクトルについてlistDistを1回だけ計算し、例えばグローバル環境内に結果リストを格納する方が良いかもしれない。

---

かいつまん

このような問題を治療するために、通常の方法では、シンプルで早いのが特長です：

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

は、私は、デフォルトでは、この方法を使用することをお勧め。メモリ限界に達する状況、すなわち非常に大きなベクトルの場合には、より複雑な解決策が必要になることがある。vec1。上述のリストベースのアプローチは、その後、救済を提供することができる。