1.Все расстояния
- только с использованием
numpy
Наивный метод:
D = np.sqrt(np.sum((X[:, None, :] - Y[None, :, :])**2, axis = -1))
Однако это занимаетмного памяти создает (i, j, n) -образную промежуточную матрицу и очень медленно
Однако, благодаря трюку от @Divakar (пакет eucl_dist, wiki ), мы можем использовать немного алгебры и np.einsum для разложения таким образом: (X - Y)**2 = X**2 - 2*X*Y + Y**2
D = np.sqrt( # (X - Y) ** 2
np.einsum('ij, ij ->i', X, X)[:, None] + # = X ** 2 \
np.einsum('ij, ij ->i', Y, Y) - # + Y ** 2 \
2 * X.dot(Y.T)) # - 2 * X * Y
Аналогично предыдущему:
XX = np.einsum('ij, ij ->i', X, X)
D = np.sqrt(XX[:, None] + XX - 2 * X.dot(X.T))
Остерегайтесь того, что неточность с плавающей запятой может привести к тому, что диагональные слагаемые будут очень незначительно отклоняться от нуля с помощью этого метода.Если вам нужно убедиться, что они равны нулю, вам нужно явно установить его:
np.einsum('ii->i', D)[:] = 0
scipy.spatial.distance.cdist является наиболее интуитивно понятной встроенной функцией для этого, и гораздо быстрее, чем просто numpy
from scipy.spatial.distance import cdist
D = cdist(X, Y)
cdist также может работать со многими, многими мерами расстояния, а также с пользователем.определенные меры расстояния (хотя они не оптимизированы).Для получения подробной информации см. Документацию, указанную выше.
Для самостоятельных расстояний scipy.spatial.distance.pdist работает аналогично cdist, но возвращает 1-D массив сжатых расстояний, экономя место на матрице симметричного расстояния, используя каждый член только один раз.Вы можете преобразовать это в квадратную матрицу, используя squareform
from scipy.spatial.distance import pdist, squareform
D_cond = pdist(X)
D = squareform(D_cond)
2.K Ближайшие соседи (KNN)
- Только с использованием
numpy
Мы могли бы использовать np.argpartition, чтобы получить индексы k-nearest и использоватьте, чтобы получить соответствующие значения расстояния.Таким образом, с D в качестве массива, содержащего значения расстояния, полученные выше, мы получили бы -
if k == 1:
k_i = D.argmin(0)
else:
k_i = D.argpartition(k, axis = 0)[:k]
k_d = np.take_along_axis(D, k_i, axis = 0)
Однако мы можем немного ускорить это, не принимая квадратные корни, пока мы не уменьшим наш набор данных.np.sqrt - самая медленная часть расчета евклидовой нормы, поэтому мы не хотим этого делать до конца.
D_sq = np.einsum('ij, ij ->i', X, X)[:, None] +\
np.einsum('ij, ij ->i', Y, Y) - 2 * X.dot(Y.T)
if k == 1:
k_i = D_sq.argmin(0)
else:
k_i = D_sq.argpartition(k, axis = 0)[:k]
k_d = np.sqrt(np.take_along_axis(D_sq, k_i, axis = 0))
Теперь np.argpartition выполняет косвенное разбиение и не обязательно дает намэлементы в отсортированном порядке и только гарантирует, что первые k элементы самые маленькие.Итак, для отсортированного вывода нам нужно использовать argsort на выходе предыдущего шага -
sorted_idx = k_d.argsort(axis = 0)
k_i_sorted = np.take_along_axis(k_i, sorted_idx, axis = 0)
k_d_sorted = np.take_along_axis(dist, k_i_sorted, axis = 0)
KD-Tree - гораздо более быстрый метод поиска соседей и ограниченных расстояний.Наиболее рекомендуемый метод - использовать scipy scipy.spatial.KDTree или scipy.spatial.cKDTree
from scipy.spatial.distance import KDTree
X_tree = KDTree(X)
k_d, k_i = X_tree.query(Y, k = k)
К сожалению, реализация scipy в KDTree идет медленнои имеет тенденцию к segfault для больших наборов данных.Как указывает @HansMusgrave здесь , pykdtree значительно повышает производительность, но не так часто, как include, как scipy и в настоящее время может работать только с евклидовым расстоянием (в то время какKDTree in scipy может обрабатывать p-нормы Минковси любого порядка)
BallTree имеет аналогичные алгоритмические свойствак KDTree.Я не знаю о параллельном / векторизованном / быстром BallTree в Python, но используя scipy, у нас все еще могут быть разумные запросы KNN для пользовательских метрик.Если доступно, встроенные метрики будут намного быстрее.
def d(a, b):
return max(np.abs(a-b))
tree = sklearn.neighbors.BallTree(X, metric=d)
k_d, k_i = tree.query(Y)
Этот ответ будет неправильным , если d() не является метрикой .Единственная причина, по которой BallTree работает быстрее, чем грубая сила, заключается в том, что свойства метрики позволяют исключить некоторые решения.Для действительно произвольных функций фактически необходима грубая сила.
3.Радиус поиска
- Только с использованием
numpy
Самый простой способ - просто использовать логическое индексирование:
mask = D_sq < r**2
r_i, r_j = np.where(mask)
r_d = np.sqrt(D_sq[mask])
Как и выше, вы можете использовать scipy.spatial.KDTree.query_ball_point
r_ij = X_tree.query_ball_point(Y, r = r)
или scipy.spatial.KDTree.query_ball_tree
Y_tree = KDTree(Y)
r_ij = X_tree.query_ball_tree(Y_tree, r = r)
К сожалению r_ij в итоге представляет собой список индексных массивов, которые немного сложно распутать для последующего использования.
Гораздо проще использовать cKDTree sparse_distance_matrix, который может выдавать coo_matrix
from scipy.spatial.distance import cKDTree
X_cTree = cKDTree(X)
Y_cTree = cKDTree(Y)
D_coo = X_cTree.sparse_distance_matrix(Y_cTree, r = r, output_type = `coo_matrix`)
r_i = D_coo.row
r_j = D_coo.column
r_d = D_coo.data
Это необычайно гибкий формат для матрицы расстояний, поскольку он остается фактической матрицей, которую (при преобразовании в csr) можно также использовать для многих векторизованных операций.