R - Как получить индексы строк и столбцов соответствующих элементов из матрицы расстояний

У меня есть целочисленный вектор 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 возражать против матрицы или перебирать матрицу?

Источник

2 ответа

Решение

Матрица расстояний - это нижняя треугольная матрица в упакованном формате, где нижняя треугольная ячейка хранится как одномерный вектор за столбцом. Вы можете проверить это через

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

Даже если мы позвоним dist(vec1, diag = TRUE, upper = TRUE) вектор все тот же; меняются только стили печати. То есть как ни звони dist Вы всегда получаете вектор.

Этот ответ сфокусирован на том, как преобразовать между 1D и 2D индексом, чтобы вы могли работать с "dist" объектом, предварительно не превратив его в полную матрицу, используя as.matrix, Если вы хотите сделать матрицу, используйте dist2mat функция, определенная в as.matrix для объекта расстояния, очень медленная; как сделать это быстрее?,

2D в 1D

1D в 2D

R функции

Для этих индексных преобразований легко написать векторизованные R-функции. Нам нужна лишь некоторая осторожность при работе с индексом "вне границ", для которого 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

Вообще говоря, матрица расстояний содержит числа с плавающей точкой. Это рискованно использовать == судить, равны ли два числа с плавающей запятой. Читать Почему эти цифры не равны? для больше и возможных стратегий.

Альтернатива с dist2mat

С использованием dist2mat функция, заданная в as.matrix для объекта расстояния, очень медленная; как сделать это быстрее? мы можем использовать 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 и использовать which с возможностью arr.ind=TRUE, Единственным недостатком этого стандартного метода для извлечения индексных чисел в матрице 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)))
  # http://stackru.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 x N. Хотя эти преобразования не тривиальны, они просты.

В комментарии в коде я указал на ответ на Stackru, который был использован здесь, чтобы разбить вектор на список. Циклы (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, Описанный выше подход, основанный на списках, мог бы тогда помочь.

Источник
Другие вопросы по тегам