我正在寻找有关快速、数值稳定的成对欧几里德距离算法的资源。特别是，假设$A \in \mathbb{R}^{M \times D}$和$B \in \mathbb{R}^{N \times D}$是两组行向量。我想计算矩阵，

$$X \in \mathbb{R}^{M \times N}, \quad X_{i,j} = \|A_i - B_j\|_2^2.$$

到目前为止，我找到了两种方法：

方法 1 - 简单的方法是遍历每个向量$A_i$和每个向量$B_j$并用相应的平方和填充$X$

方法 2 - 利用以下事实，

$$\|A_i - B_j\|_2^2 = \langle A_i - B_j, A_i - B_j \rangle = \|A_i\|_2^2 + \|B_j\|_2^2 - 2 \langle A_i, B_j \rangle.$$

利用高效的代码来计算这三个项中的每一项，方法 2 几乎快了一个数量级。但是，它在数值上不如方法 1 稳定。例如，方法 2 可以输出负距离。

是否有替代方法来计算比方法 1 更快但保证（至少）所有非负距离的成对距离？有更好的数值稳定性？

代码演示- 下面是一个用 Python 编写的示例比较。Scipy 的cdist()功能实际上是方法 1 的实现，而cdist_fast()下面是方法 2 的实现：

# experiment.py
import numpy as np
import time
from scipy.spatial.distance import cdist


def cdist_fast(XA, XB):
    XA_norm = np.sum(XA**2, axis=1)
    XB_norm = np.sum(XB**2, axis=1)
    XA_XB_T = np.dot(XA, XB.T)
    distances = XA_norm.reshape(-1,1) + XB_norm - 2*XA_XB_T
    return distances


def main():
    M,N = 5000, 128
    XA = np.random.randn(M,N)

    t = time.time()
    distances_cdist = cdist(XA, XA, metric='sqeuclidean')
    time_cdist = time.time() - t

    t = time.time()
    distances_cdist_fast = cdist_fast(XA, XA)
    time_cdist_fast = time.time() - t

    print(f'time_cdist = {time_cdist:.3f} s')
    print(f'time_cdist_fast = {time_cdist_fast:.3f} s')

    # check validity of results
    assert np.allclose(distances_cdist, distances_cdist_fast)

    # check that the results are non-negative
    try:
        assert (distances_cdist >= 0.0).all()
    except AssertionError:
        print('Numerical instability in cdist()')

    try:
        assert (distances_cdist_fast >= 0.0).all()
    except AssertionError:
        print('Numerical instability in cdist_fast()')


if __name__ == '__main__':
    main()

3.1 GHz Intel Core i7 上的脚本输出：

$ python experiment.py
time_cdist = 3.457 s
time_cdist_fast = 0.625 s
Numerical instability in cdist_fast()