How can the Euclidean distance be calculated with NumPy?

У меня есть две точки в трехмерном пространстве:

a = (ax, ay, az)
b = (bx, by, bz)

Я хочу вычислить расстояние между ними:

dist = sqrt((ax-bx)^2 + (ay-by)^2 + (az-bz)^2)

Как мне это сделать с помощью NumPy? У меня есть:

import numpy
a = numpy.array((ax, ay, az))
b = numpy.array((bx, by, bz))

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

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

Это работает, потому что евклидово расстояние является нормой l2, а значение по умолчанию ord параметра в numpy.linalg.norm равно 2.
Подробнее о теории читайте в Введении в интеллектуальный анализ данных:

Использование scipy.spatial.distance.euclidean:

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

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

Первый совет - организовать ваши данные таким образом, чтобы массивы имели размерность (3, n) (и, очевидно, были C-смежными). Если добавление происходит в непрерывном первом измерении, все происходит быстрее, и не имеет особого значения, используете ли вы sqrt-sum with axis=0, linalg.norm with axis=0 или

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

это, с небольшим отрывом, самый быстрый вариант. (На самом деле это справедливо и для одной строки.)

Варианты, в которых вы суммируете по второй оси, axis=1, все существенно медленнее.

Код для воспроизведения графика:

import numpy
import perfplot
from scipy.spatial import distance


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


def linalg_norm_T(data):
    a, b = data[1]
    return numpy.linalg.norm(a - b, axis=0)


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


def sqrt_sum_T(data):
    a, b = data[1]
    return numpy.sqrt(numpy.sum((a - b) ** 2, axis=0))


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


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


def sqrt_einsum_T(data):
    a, b = data[1]
    a_min_b = a - b
    return numpy.sqrt(numpy.einsum("ij,ij->j", a_min_b, a_min_b))


def setup(n):
    a = numpy.random.rand(n, 3)
    b = numpy.random.rand(n, 3)
    out0 = numpy.array([a, b])
    out1 = numpy.array([a.T, b.T])
    return out0, out1


b = perfplot.bench(
    setup=setup,
    n_range=[2 ** k for k in range(22)],
    kernels=[
        linalg_norm,
        linalg_norm_T,
        scipy_distance,
        sqrt_sum,
        sqrt_sum_T,
        sqrt_einsum,
        sqrt_einsum_T,
    ],
    xlabel="len(x), len(y)",
)
b.save("norm.png")

Я хочу разъяснить простой ответ с различными замечаниями по производительности. 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 не кэширует поиск по именам.

Итак

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)

Первое, что нам нужно помнить, это то, что мы используем Pythagoras для вычисления расстояния (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)

Отлично, обе функции больше не извлекают дорогостоящих квадратных корней. Это будет намного быстрее, но прежде чем идти дальше, проверьте себя: зачем sort_things_by_distance понадобился "наивный" отказ от ответственности оба раза выше? Ответ в самом низу (* a1).

Мы можем улучшить 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))

(* a1: ключ сортировки sort_things_by_distance вызывает distance_sq для каждого отдельного элемента, и этот невинно выглядящий ключ является лямбдой , которая является второй функцией, которая должна быть вызвана ...)