scaled/rotated 成对平方欧氏距离的矢量化计算

Question

给定一组 n 维度的向量 d 存储在 (n,d) 数组和第二组相同维度的 m 向量（存储在 (m,d) 数组中）我想计算平方向量之间的逐点距离，由大小为 (d,d).

的矩阵 A 缩放

输出应该是一个(n,m)数组。

我希望 m 和 n 的输入范围在 1 到 10.000 之间，[=38= 的输入范围在 1 到 100 之间]d.

两点之间的距离由下式给出：

在未优化但有效的 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 中的嵌套循环，例如使用广播黑魔法？

编辑：

案例

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

这是正常的平方欧几里德距离，可以计算例如

from scipy.spatial.distance import cdist

cdist(v1,v2,'sqeuclidean')

Answer 1

我们可以使用np.einsum-

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

此外，通过将 np.einsum 中的 optimize 标志设置为 True 来使用 BLAS。

矢量化方法说明

原代码是-

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

我。我们正在翻译：

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

有关 broadcasting 的更多信息可在 docs 中找到。

二．接下来，(v1[i,:] - v2[j,:]).T @ A @ (v1[i,:] - v2[j,:]) 和新的 V 使用 einsum 的字符串表示法变为 np.einsum('ijk,kl,ijl->ij',V,A,V)。更多信息可以在 docs.

中找到

scaled/rotated 成对平方欧氏距离的矢量化计算

Vectorized calculation of scaled/rotated pairwise squared euclidean distance

python

numpy

vectorization

array-broadcasting