Как наиболее эффективно рассчитать квадратное евклидово расстояние между N образцами и центроидами кластеров?

Question

Я ищу эффективный способ ( нет для петель ) для вычисления евклидова расстояния между набором выборок и набором центроидов кластеров.

Пример:

import numpy as np
X = np.array([[1,2,3],[1, 1, 1],[0, 2, 0]])
y = np.array([[1,2,3], [0, 1, 0]])

Ожидаемый результат:

array([[ 0., 11.],
       [ 5.,  2.],
       [10.,  1.]])

Это квадрат евклидова расстояния между каждым образцом в X и каждым центроидом в y.

Я придумал 2 решения:

Решение 1:

def dist_2(X,y):
    X_square_sum = np.sum(np.square(X), axis = 1)
    y_square_sum = np.sum(np.square(y), axis = 1)
    dot_xy = np.dot(X, y.T)
    X_square_sum_tile = np.tile(X_square_sum.reshape(-1, 1), (1, y.shape[0]))
    y_square_sum_tile = np.tile(y_square_sum.reshape(1, -1), (X.shape[0], 1))
    dist = X_square_sum_tile + y_square_sum_tile - (2 * dot_xy)
    return dist

dist = dist_2(X, y)

решение 2:

import scipy
dist = scipy.spatial.distance.cdist(X,y)**2

Производительность (настенное время) решения «два»

import time
X = np.random.random((100000, 50))
y = np.random.random((100, 50))

start = time.time()
dist = scipy.spatial.distance.cdist(X,y)**2
end = time.time()
print (end - start)

Среднее время нахождения настенных часов = 0,7 с

start = time.time()
dist = dist_2(X,y)
end = time.time()
print (end - start)

Среднее время нахождения настенных часов = 0,3 с

Тест на большое количество центроидов

X = np.random.random((100000, 50))
y = np.random.random((1000, 50))

Среднее время нахождения настенных часов в «решении 1» = 50 с (+ проблема с памятью)

Среднее время нахождения настенных часов в «решении 2» = 6 секунд !!!

Заключение

Похоже, что "решение 1 более эффективно, чем" решение 2 "в отношении среднего истекшего времени настенных часов (для небольших наборов данных), но неэффективно в отношении памяти.

Есть предложения?

Answer 1

Этот вопрос часто задают в сочетании с поиском ближайших соседей.Если это так, взгляните на подход kdtree .Это будет гораздо более эффективно, чем вычисление евклидовых расстояний, как по потреблению памяти, так и по производительности.

Если это не так, вот несколько возможных подходов.Первые два из являются ответом Divakar .Третий использует Numba для JIT-компиляции.Основное отличие первых двух версий состоит в том, что он избегает временных массивов.

Три подхода к вычислению евклидовых расстояний

import numpy as np
import numba as nb

# @Paul Panzer
#https://stackoverflow.com/a/42994680/4045774
def outer_sum_dot_app(A,B):
    return np.add.outer((A*A).sum(axis=-1), (B*B).sum(axis=-1)) - 2*np.dot(A,B.T)

# @Divakar
#https://stackoverflow.com/a/42994680/4045774
def outer_einsum_dot_app(A,B):
    return np.einsum('ij,ij->i',A,A)[:,None] + np.einsum('ij,ij->i',B,B) - 2*np.dot(A,B.T)

@nb.njit()
def calc_dist(A,B,sqrt=False):
  dist=np.dot(A,B.T)

  TMP_A=np.empty(A.shape[0],dtype=A.dtype)
  for i in range(A.shape[0]):
    sum=0.
    for j in range(A.shape[1]):
      sum+=A[i,j]**2
    TMP_A[i]=sum

  TMP_B=np.empty(B.shape[0],dtype=A.dtype)
  for i in range(B.shape[0]):
    sum=0.
    for j in range(B.shape[1]):
      sum+=B[i,j]**2
    TMP_B[i]=sum

  if sqrt==True:
    for i in range(A.shape[0]):
      for j in range(B.shape[0]):
        dist[i,j]=np.sqrt(-2.*dist[i,j]+TMP_A[i]+TMP_B[j])
  else:
    for i in range(A.shape[0]):
      for j in range(B.shape[0]):
        dist[i,j]=-2.*dist[i,j]+TMP_A[i]+TMP_B[j]
  return dist

Время

A = np.random.randn(10000,3)
B = np.random.randn(10000,3)

#calc_dist:                      360ms first call excluded due to compilation overhead
#outer_einsum_dot_app (Divakar): 1150ms
#outer_sum_dot_app (Paul Panzer):1590ms
#dist_2:                         1840ms

A = np.random.randn(1000,100)
B = np.random.randn(1000,100)

#calc_dist:                      4.3  ms first call excluded due to compilation overhead
#outer_einsum_dot_app (Divakar): 13.12ms
#outer_sum_dot_app (Paul Panzer):13.2 ms
#dist_2:                         21.3ms