Как рассчитать евклидово расстояние с помощью NumPy?

Question

У меня есть две точки в 3D:

(xa, ya, za)
(xb, yb, zb)

И я хочу рассчитать расстояние:

dist = sqrt((xa-xb)^2 + (ya-yb)^2 + (za-zb)^2)

Какой лучший способ сделать это с NumPy или с Python в целом? У меня есть:

a = numpy.array((xa ,ya, za))
b = numpy.array((xb, yb, zb))

Answer 1

Использование numpy.linalg.norm:

dist = numpy.linalg.norm(a-b)

Answer 2

В SciPy есть функция для этого. Это называется Евклидов .

Пример:

from scipy.spatial import distance
a = (1, 2, 3)
b = (4, 5, 6)
dst = distance.euclidean(a, b)

Answer 3

Для тех, кто заинтересован в вычислении нескольких расстояний одновременно, я провел небольшое сравнение, используя perfplot (мой маленький проект). Оказывается,

a_min_b = a - b
numpy.sqrt(numpy.einsum('ij,ij->i', a_min_b, a_min_b))

вычисляет расстояния строк в a и b быстрее всего. Это справедливо и для одной строки!

---

Код для воспроизведения сюжета:

import matplotlib
import numpy
import perfplot
from scipy.spatial import distance


def linalg_norm(data):
    a, b = data
    return numpy.linalg.norm(a-b, axis=1)


def sqrt_sum(data):
    a, b = data
    return numpy.sqrt(numpy.sum((a-b)**2, axis=1))


def scipy_distance(data):
    a, b = data
    return list(map(distance.euclidean, a, b))


def mpl_dist(data):
    a, b = data
    return list(map(matplotlib.mlab.dist, a, b))


def sqrt_einsum(data):
    a, b = data
    a_min_b = a - b
    return numpy.sqrt(numpy.einsum('ij,ij->i', a_min_b, a_min_b))


perfplot.show(
    setup=lambda n: numpy.random.rand(2, n, 3),
    n_range=[2**k for k in range(20)],
    kernels=[linalg_norm, scipy_distance, mpl_dist, sqrt_sum, sqrt_einsum],
    logx=True,
    logy=True,
    xlabel='len(x), len(y)'
    )

Answer 4

Другой экземпляр этого метода решения проблем :

def dist(x,y):   
    return numpy.sqrt(numpy.sum((x-y)**2))

a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))
dist_a_b = dist(a,b)

Answer 5

Я хочу изложить простой ответ с различными заметками о производительности. np.linalg.norm сделает, возможно, больше, чем нужно:

dist = numpy.linalg.norm(a-b)

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

sP = set(points)
pA = point
distances = np.linalg.norm(sP - pA, ord=2, axis=1.)  # 'distances' is a list

Помните несколько вещей:

• Вызовы функций Python стоят дорого.
• [Обычный] Python не кэширует поиск имен.

So

def distance(pointA, pointB):
    dist = np.linalg.norm(pointA - pointB)
    return dist

не так невинен, как выглядит.

>>> dis.dis(distance)
  2           0 LOAD_GLOBAL              0 (np)
              2 LOAD_ATTR                1 (linalg)
              4 LOAD_ATTR                2 (norm)
              6 LOAD_FAST                0 (pointA)
              8 LOAD_FAST                1 (pointB)
             10 BINARY_SUBTRACT
             12 CALL_FUNCTION            1
             14 STORE_FAST               2 (dist)

  3          16 LOAD_FAST                2 (dist)
             18 RETURN_VALUE

Во-первых - каждый раз, когда мы вызываем его, мы должны выполнить глобальный поиск для «np», поиск в области для «linalg» и поиск в области для «norm», а также накладные расходы, связанные с вызовом функция может приравниваться к десяткам инструкций Python.

Наконец, мы потратили две операции, чтобы сохранить результат и перезагрузить его для возврата ...

Первый проход при улучшении: сделайте поиск быстрее, пропустите магазин

def distance(pointA, pointB, _norm=np.linalg.norm):
    return _norm(pointA - pointB)

Мы получаем гораздо более обтекаемый:

>>> dis.dis(distance)
  2           0 LOAD_FAST                2 (_norm)
              2 LOAD_FAST                0 (pointA)
              4 LOAD_FAST                1 (pointB)
              6 BINARY_SUBTRACT
              8 CALL_FUNCTION            1
             10 RETURN_VALUE

Затраты на вызов функции по-прежнему составляют некоторую работу. И вы захотите сделать тесты, чтобы определить, лучше ли вам делать математику самостоятельно:

def distance(pointA, pointB):
    return (
        ((pointA.x - pointB.x) ** 2) +
        ((pointA.y - pointB.y) ** 2) +
        ((pointA.z - pointB.z) ** 2)
    ) ** 0.5  # fast sqrt

На некоторых платформах **0.5 быстрее, чем math.sqrt. Ваш пробег может отличаться.

**** Расширенные замечания по производительности.

Почему вы рассчитываете расстояние? Если единственной целью является его отображение,

 print("The target is %.2fm away" % (distance(a, b)))

двигаться вперед. Но если вы сравниваете расстояния, проводите проверки дальности и т. Д., Я хотел бы добавить некоторые полезные наблюдения за производительностью.

Давайте рассмотрим два случая: сортировка по расстоянию или отбор списка по элементам, которые соответствуют ограничению диапазона.

# Ultra naive implementations. Hold onto your hat.

def sort_things_by_distance(origin, things):
    return things.sort(key=lambda thing: distance(origin, thing))

def in_range(origin, range, things):
    things_in_range = []
    for thing in things:
        if distance(origin, thing) <= range:
            things_in_range.append(thing)

Первое, что нам нужно запомнить, это то, что мы используем Пифагор для вычисления расстояния (dist = sqrt(x^2 + y^2 + z^2)), поэтому мы делаем много sqrt вызовов. Математика 101:

dist = root ( x^2 + y^2 + z^2 )
:.
dist^2 = x^2 + y^2 + z^2
and
sq(N) < sq(M) iff M > N
and
sq(N) > sq(M) iff N > M
and
sq(N) = sq(M) iff N == M

Короче говоря: до тех пор, пока нам не потребуется расстояние в единице X, а не X ^ 2, мы можем исключить самую сложную часть вычислений.

# Still naive, but much faster.

def distance_sq(left, right):
    """ Returns the square of the distance between left and right. """
    return (
        ((left.x - right.x) ** 2) +
        ((left.y - right.y) ** 2) +
        ((left.z - right.z) ** 2)
    )

def sort_things_by_distance(origin, things):
    return things.sort(key=lambda thing: distance_sq(origin, thing))

def in_range(origin, range, things):
    things_in_range = []

    # Remember that sqrt(N)**2 == N, so if we square
    # range, we don't need to root the distances.
    range_sq = range**2

    for thing in things:
        if distance_sq(origin, thing) <= range_sq:
            things_in_range.append(thing)

Отлично, обе функции больше не делают дорогих квадратных корней. Это будет намного быстрее. Мы также можем улучшить in_range, преобразовав его в генератор:

def in_range(origin, range, things):
    range_sq = range**2
    yield from (thing for thing in things
                if distance_sq(origin, thing) <= range_sq)

Это особенно полезно, если вы делаете что-то вроде:

if any(in_range(origin, max_dist, things)):
    ...

Но если следующая вещь, которую вы собираетесь сделать, требует расстояния,

for nearby in in_range(origin, walking_distance, hotdog_stands):
    print("%s %.2fm" % (nearby.name, distance(origin, nearby)))

рассмотрите возможность получения кортежей:

def in_range_with_dist_sq(origin, range, things):
    range_sq = range**2
    for thing in things:
        dist_sq = distance_sq(origin, thing)
        if dist_sq <= range_sq: yield (thing, dist_sq)

Это может быть особенно полезно, если вы можете чередовать проверки диапазона («находите вещи, которые находятся около X и в пределах Nm от Y», поскольку вам не нужно снова вычислять расстояние).

Но что если мы будем искать действительно большой список things, и мы ожидаем, что многие из них не заслуживают рассмотрения?

На самом деле очень простая оптимизация:

def in_range_all_the_things(origin, range, things):
    range_sq = range**2
    for thing in things:
        dist_sq = (origin.x - thing.x) ** 2
        if dist_sq <= range_sq:
            dist_sq += (origin.y - thing.y) ** 2
            if dist_sq <= range_sq:
                dist_sq += (origin.z - thing.z) ** 2
                if dist_sq <= range_sq:
                    yield thing

Будет ли это полезно, зависит от размера "вещей".

def in_range_all_the_things(origin, range, things):
    range_sq = range**2
    if len(things) >= 4096:
        for thing in things:
            dist_sq = (origin.x - thing.x) ** 2
            if dist_sq <= range_sq:
                dist_sq += (origin.y - thing.y) ** 2
                if dist_sq <= range_sq:
                    dist_sq += (origin.z - thing.z) ** 2
                    if dist_sq <= range_sq:
                        yield thing
    elif len(things) > 32:
        for things in things:
            dist_sq = (origin.x - thing.x) ** 2
            if dist_sq <= range_sq:
                dist_sq += (origin.y - thing.y) ** 2 + (origin.z - thing.z) ** 2
                if dist_sq <= range_sq:
                    yield thing
    else:
        ... just calculate distance and range-check it ...

И снова, рассмотрите возможность выдачи dist_sq. Наш пример с хот-догом становится:

# Chaining generators
info = in_range_with_dist_sq(origin, walking_distance, hotdog_stands)
info = (stand, dist_sq**0.5 for stand, dist_sq in info)
for stand, dist in info:
    print("%s %.2fm" % (stand, dist))

Answer 6

Я нахожу функцию dist в matplotlib.mlab, но я не думаю, что она достаточно удобна.

Я публикую это здесь только для справки.

import numpy as np
import matplotlib as plt

a = np.array([1, 2, 3])
b = np.array([2, 3, 4])

# Distance between a and b
dis = plt.mlab.dist(a, b)

Answer 7

Это можно сделать следующим образом. Я не знаю, как это быстро, но он не использует NumPy.

from math import sqrt
a = (1, 2, 3) # Data point 1
b = (4, 5, 6) # Data point 2
print sqrt(sum( (a - b)**2 for a, b in zip(a, b)))

Answer 8

Мне нравится np.dot (точечный продукт):

a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))

distance = (np.dot(a-b,a-b))**.5

Answer 9

Вы можете просто вычесть векторы и затем произвести их.

Следуя вашему примеру,

a = numpy.array((xa, ya, za))
b = numpy.array((xb, yb, zb))

tmp = a - b
sum_squared = numpy.dot(tmp.T, tmp)
result sqrt(sum_squared)

Это простой код и его легко понять.

Answer 10

Начиная с Python 3.8, модуль math напрямую предоставляет функцию dist, которая возвращает евклидово расстояние между двумя точками (в виде набора координат):

from math import dist

dist((1, 2, 6), (-2, 3, 2)) # 5.0990195135927845

Если вы работаете со списками вместо кортежей:

dist(tuple([1, 2, 6]), tuple([-2, 3, 2]))