Вычислить евклидово расстояние на numpy перекрестном произведении строка-строка?

Question

У меня есть два массива numpy с n количеством координат (два элемента в строке).

coords_a = np.random.random((20, 2))
coords_b = np.random.random((20, 2))

Теперь для каждой комбинации строк я хочу вычислить функцию и сохранить возвращаемое значение как элемент в матрице. Таким образом, результирующий массив должен иметь форму (20, 20) и его можно "лениво" вычислять, как показано ниже. В качестве примерной функции используется евклидово расстояние.

def euclidean_dist(x1: float, y1: float, x2: float, y2: float) -> float:
    """Return the euclidean distance between two the points (x1, y1) and (x2, y2)."""
    return np.sqrt(np.square(x1 - x2) + np.square(y1 - y2))

matrix = []
for a in coords_a:
    row = []
    for b in coords_b:
        row.append(euclidean_dist(*a, *b))
    matrix.append(row)
    
matrix = np.array(matrix)

Как вы можете себе представить, эта вложенность для l oop занимает очень много времени, занимая более 25 секунд всего с 2000 пар координат. Есть ли рекомендуемый способ векторизации такого перекрестного произведения?

Заранее спасибо.

Answer 1

Я хотел бы добавить свои 2 цента, поскольку не все функции уже реализованы в numpy или scipy. В общем, вы можете использовать numpy широковещательную передачу для получения векторизованного решения. Для конкретного c случая евклидова расстояния здесь, как вы это делаете:

import numpy as np

# Define the arrays of coordinates
coords_a = np.random.random((20, 2))
coords_b = np.random.random((20, 2))

# Expand their dimensions
a = coords_a[:, None]
b = coords_b[None, None]

# Use broadcasting to compute pairwise difference
d = a-b

# Apply formula for euclidean distance
r = np.sqrt(np.sum(d**2, axis=-1)) 

С точки зрения производительности времени для этого c случай scipy.spatial.distance.cdist намного быстрее, но не все функции доступны :

import numpy as np
from scipy.spatial.distance import cdist

a = np.random.random((10_000, 2))
b = np.random.random((10_000, 2))

euc_broadcast = lambda a,b: np.sqrt(np.sum(np.square(a[:, None]-b[None, :]), axis=-1))

%timeit euc_broadcast(a, b)
3.39 s ± 149 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit cdist(a, b)
603 ms ± 13.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Answer 2

Для вашего конкретного примера c вы можете сделать:

from scipy.spatial.distance import cdist
cdist(coords_b,coords_a)

В общем, векторизация зависит от вашей функции.