R – Come ottenere gli indici di riga e colonna di elementi abbinati da una matrice di distanze

Ho un vettore intero vec1 e sto generando una matrice distante usando la funzione dist . Voglio ottenere le coordinate (riga e colonna) dell’elemento di un determinato valore nella matrice della distanza. Essenzialmente vorrei ottenere la coppia di elementi distanti tra loro. Per esempio:

 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 

Dite, mi interessa la coppia di elementi nel vettore che sono 5 unità a parte. Volevo ottenere la coordinata1 che sono le righe e le coordinate2 che sono le colonne della matrice della distanza. In questo esempio di giocattolo, mi aspetterei

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

Mi chiedo se esiste un modo efficace per ottenere questi valori che non implichino la conversione dell’object dist in una matrice o il looping della matrice?

Una matrice di distanze è una matrice triangular inferiore in formato compatto, in cui il triangular inferiore viene memorizzato come vettore 1D per colonna. Puoi controllare questo via

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

Anche se chiamiamo dist(vec1, diag = TRUE, upper = TRUE) , il vettore è sempre lo stesso; cambiano solo gli stili di stampa. Cioè, non importa come chiami dist , ottieni sempre un vettore.

Questa risposta si concentra su come trasformare tra l’indice 1D e 2D, in modo da poter lavorare con un object “dist” senza prima trasformarlo in una matrice completa usando as.matrix . Se si desidera renderlo una matrice, utilizzare la funzione dist2mat definita in as.matrix su un object distanza è estremamente lenta; come renderlo più veloce? .


Da 2D a 1D

Da 1D a 2D


Funzioni R

È facile scrivere funzioni R vettorizzate per quelle trasformazioni di indice. Abbiamo solo bisogno di cure per l’indice “out-of-bound”, per il quale NA deve essere restituito.

 ## 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) } 

Queste funzioni sono estremamente economiche nell'utilizzo della memoria, poiché funzionano con indici anziché matrici.


Applicando finv alla tua domanda

Puoi usare

 vec1 <- c(2,3,6,12,17) distMatrix <- dist(vec1) finv(which(distMatrix == 5), distMatrix) # ij #[1,] 5 4 

In generale, una matrice di distanze contiene numeri in virgola mobile. È rischioso usare == per giudicare se due numeri in virgola mobile sono uguali. Leggi Perché questi numeri non sono uguali? per più e possibili strategie.


Alternativa con dist2mat

L'uso della funzione dist2mat fornita in as.matrix su un object distanza è estremamente lenta; come renderlo più veloce? , possiamo usare 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 

Appendice: Markdown (richiede supporto MathJax) per l'immagine

 ## 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 jth 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)$th "$\times$" in the $j$th 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$th "$\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$th column, ie, $(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 (ie, 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\}.$$ 

Se il vettore non è troppo grande, il modo migliore è probabilmente di avvolgere l’output di dist in as.matrix e di usare which con l’opzione arr.ind=TRUE . L’unico svantaggio di questo metodo standard per recuperare i numeri indice all’interno di una dist distinta è un aumento dell’uso della memoria, che può diventare importante nel caso di vettori molto grandi passati a dist . Questo perché la conversione della matrice triangular inferiore restituita da dist in una matrice densa e regolare raddoppia efficacemente la quantità di dati memorizzati.

Un’alternativa consiste nel convertire l’object dist in una lista, in modo tale che ciascuna colonna nella matrice triangular inferiore di dist rappresenti un membro della lista. Il numero di indice dei membri dell’elenco e la posizione degli elementi all’interno dei membri dell’elenco possono quindi essere associati alla colonna e al numero di riga della matrice N x N densa, senza generare la matrice.

Ecco una ansible implementazione di questo approccio basato su elenchi:

 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 

Le parti del codice che potrebbero essere alquanto poco chiare riguardano la mapping della posizione di una voce all'interno della lista a un valore di colonna / riga della matrice N x N. Anche se non banali, queste trasformazioni sono semplici.

In un commento all'interno del codice ho indicato una risposta su StackOverflow che è stata usata qui per dividere un vettore in una lista. I loop (sapply, lapply) non dovrebbero essere problematici in termini di prestazioni poiché il loro range è dell'ordine O (N). L'utilizzo della memoria di questo codice è in gran parte determinato dalla memorizzazione della lista. Questa quantità di memoria dovrebbe essere simile a quella dell'object dist poiché entrambi gli oggetti contengono gli stessi dati.

L'object dist viene calcolato e trasformato in una lista nella funzione distToList() . A causa del calcolo dist, che è richiesto in ogni caso, questa funzione potrebbe richiedere molto tempo nel caso di vettori di grandi dimensioni. Se l'objective è trovare più coppie con valori di distanza diversi, allora potrebbe essere meglio calcolare listDist solo una volta per un dato vettore e memorizzare l'elenco risultante, ad esempio, nell'ambiente globale.


Per farla breve

Il modo usuale per trattare tali problemi è semplice e veloce:

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

Suggerisco di utilizzare questo metodo per impostazione predefinita. Soluzioni più complicate possono rendersi necessarie in situazioni in cui si raggiungono limiti di memoria, cioè, nel caso di vettori molto grandi, vec1 . L'approccio basato su elenchi descritto sopra potrebbe quindi fornire un rimedio.