计算（平方）欧氏距离的替代方法的数值问题

Question

我想快速计算平方欧几里德，如下所述：

注1：我只对距离感兴趣，不是RBF核

注2：我这里忽略了numexpr，直接用了numpy

简而言之，我计算：

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

我计算距离矩阵的速度比 scipy.pdist 快 ~10 倍。但是，我观察到数值问题，如果我取平方根来获得欧氏距离，这些问题会变得更糟。我的值按 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]

正如您所看到的，一些值也是负数（在采用 sqrt 后结果为 nan）。当然这些很容易捕捉——但对于欧几里得情况来说，小的积极因素有很大的错误（例如 abs_error=5.96046448e-08）

Answer 1

根据我的评论，使用 abs 可能是清除该算法固有的数值稳定性的最佳选择。由于您担心性能，因此您可能应该使用变异赋值运算符，因为它们会导致创建更少的垃圾，因此速度会更快。此外，当运行这具有许多功能（例如 10k）时，我看到 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)

在 IPython 中计时：

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

没有 errors/warnings

链接的问题还指出 scikit-learn 具有良好的距离内核良好实现，欧几里德的是 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)

如果您已经在使用该软件包，可能会更好用

计算（平方）欧氏距离的替代方法的数值问题

Numerical issues for alternative way to compute (squared) euclidean distance

python

numpy

distance

numerical-stability

问题：

示例：

示例输出：