Векторизованный расчет масштабированного / повернутого попарно квадратного евклидова расстояния - PullRequest
Новая
       33

Векторизованный расчет масштабированного / повернутого попарно квадратного евклидова расстояния

1 голос
/ 25 апреля 2019

При заданном наборе n векторов размерности d , хранящихся в массиве (n, d) и втором наборе m векторы одного и того же размера (хранятся в массиве (m, d) ). Я хочу вычислить квадратное точечное расстояние между векторами, масштабированное по некоторой матрице A с размером (d, d ) .

Выходные данные должны быть массивом (n, m) .

Я ожидаю, что диапазон ввода будет где-то между 1 и 10.000 для м и n и от 1 до 100 для d .

Расстояние между двумя точками определяется как:

d_i_j = (v1_i - v2_j)^T A (v1_i - v2_j)

В неоптимизированном, но работающем коде Python это выглядит так: 

import numpy as np

v1 = np.array([[1, 2],
               [3, 4],
               [4, 5]])

v2 = np.array([[1,1],
               [2, 2],
               [2, 2],
               [0, 0]])

A = np.array([[1,0], [2, 3]])

d = np.zeros((3, 4))

for i in range(0,3):
    for j in range(0,4):
        d[i,j] = (v1[i,:] - v2[j,:]).T @ A @ (v1[i,:] - v2[j,:])

Квадратное расстояние между точками примера: 

d = [[  3.   1.   1.  17.]
 [ 43.  17.  17.  81.]
 [ 81.  43.  43. 131.]]

Существует ли версия этого, которая позволяет избежать вложенного цикла в python, например используя трансляцию черной магии?

EDIT:

Для дела 

A = np.array([[1,0], [0, 1]])

это нормальное евклидово расстояние в квадрате, которое можно рассчитать, например: 

from scipy.spatial.distance import cdist

cdist(v1,v2,'sqeuclidean')

1 Ответ

1 голос
/ 26 апреля 2019

Мы можем использовать np.einsum - 

V = v1[:,None,:]-v2
d_out = np.einsum('ijk,kl,ijl->ij',V,A,V)

Кроме того, поиграйте с флагом optimize в np.einsum, установив для него значение True для использования BLAS.

Пояснения к векторизованному методу

Оригинальный код был - 

d[i,j] = (v1[i,:] - v2[j,:]).T @ A @ (v1[i,:] - v2[j,:])

I. Мы переводим: 

v1[i,:] - v2[j,:]

на внешнюю операцию с broadcasting: 

v1[:,None,:]-v2

Схематически обозначено: 

v1[:,None,:]  :  m x 1 x n
v2            :      m x n
output, V     :  m x m x n

Подробнее о outer объяснение.

Более подробную информацию о broadcasting можно найти в docs.

II. Затем (v1[i,:] - v2[j,:]).T @ A @ (v1[i,:] - v2[j,:]) с новым V становится np.einsum('ijk,kl,ijl->ij',V,A,V) с использованием строковой нотации einsum. Более подробную информацию можно найти в docs.

...