Численные вопросы для альтернативного способа вычисления (квадрата) евклидова расстояния

Question

Я хочу быстро вычислить квадрат евклидова, как описано здесь:

Какой самый быстрый способ вычислить ядро ​​RBF в python?

Примечание 1 : меня интересует только расстояние, не ядро ​​RBF.

Примечание2 : здесь я пренебрегаю numexpr и использовать только numpy напрямую.

Короче говоря, я вычисляю:

|| x - y ||^2 = ||x||^2 + ||y||^2 - 2. * (x @ y.T)

Я могу вычислить матрицу расстояний с коэффициентом ~10 по сравнению с scipy.pdist с этим. Однако я наблюдаю численные проблемы, которые усугубляются, если я возьму квадрат root, чтобы получить евклидово расстояние. У меня есть значения, которые имеют порядок 1E-8 - 1E-7, который должен быть ровно нулем (т.е. дублированные точки или расстояние до собственной точки).

Вопрос:

Существуют ли способы или идеи для преодоления этих числовых проблем (целесообразно, не жертвуя слишком большой скоростью оценки)? Или численные проблемы являются причиной, по которой этот путь не берется (например, scipy.pdist) в первую очередь?

Пример:

Это небольшой пример кода, показывающий числовые проблемы (не ускорения, пожалуйста, посмотрите ответы в связанной ветке SO выше).

import numpy as np

M = np.random.rand(1000, 10)

M_norm = np.sum(M**2, axis=1)

res = M_norm[:, np.newaxis] + M_norm[np.newaxis, :] - 2. * M @ M.T
unique = np.unique(np.diag(res))  # analytically all diag values are exactly zero 
sqrt_unique = np.sqrt(unique)


print(unique)
print(sqrt_unique)

Пример вывода:

[-2.66453526e-15 -1.77635684e-15 -8.88178420e-16 -4.44089210e-16
  0.00000000e+00  4.44089210e-16  8.88178420e-16  1.77635684e-15
  3.55271368e-15]
[           nan            nan            nan            nan
 0.00000000e+00 2.10734243e-08 2.98023224e-08 4.21468485e-08
 5.96046448e-08]

Как видите, некоторые значения также отрицательны (что приводит к nan после взятия sqrt). Конечно, их легко поймать - но небольшие положительные значения имеют большую ошибку для евклидова случая (например, abs_error=5.96046448e-08)

Answer 1

согласно моему комментарию, использование abs, вероятно, является лучшим вариантом для очистки числовой стабильности, присущей этому алгоритму. поскольку вы беспокоитесь о производительности, вам, вероятно, следует использовать изменяющиеся операторы присваивания, поскольку они создают меньше мусора и, следовательно, могут быть намного быстрее. Кроме того, при выполнении этого со многими функциями (например, 10 КБ) я вижу, что pdist медленнее, чем эта реализация.

, сложив вышеприведенное вместе, мы получаем:

import numpy as np

def edist0(M):
    "calculate pairwise euclidean distance"
    M_norm = np.sum(M**2, axis=1)
    res = M_norm[:, np.newaxis] + M_norm[np.newaxis, :] - 2. * M @ M.T    
    return np.sqrt(np.abs(res))

def edist1(M):
    "optimised calculation of pairwise euclidean distance"
    M_norm = np.einsum('ij,ij->i', M, M)
    res = M @ M.T
    res *= -2.
    res += M_norm[:, np.newaxis]
    res += M_norm[np.newaxis, :]
    return np.sqrt(np.abs(res, out=res), out=res)

синхронизация этого в I Python с:

from scipy.spatial import distance

M = np.random.rand(1000, 10000)

%timeit distance.squareform(distance.pdist(M))
%timeit edist0(M)
%timeit edist1(M)

Я получаю:

2.82 s ± 60.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
296 ms ± 6.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
153 ms ± 1.58 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

и без ошибок / предупреждений от sqrt

связанный вопрос также указывает на scikit-learn as с хорошими реализациями ядра расстояния, евклидовым является pairwise_distances, который оценивается как:

from sklearn.metrics import pairwise_distances

%timeit pairwise_distances(M)

170 ms ± 5.51 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

, что может быть удобно использовать, если вы уже используете этот пакет