Einsum和距离计算

Question

Einsum和距离计算

8

我正在寻找使用einsum来确定不同行数但列数相等的numpy数组之间的距离的解决方案。我尝试了各种组合，但唯一成功的方法是使用以下代码。显然我错过了什么，文献和众多线程也没有让我更接近解决方案。我希望能找到一种通用性，使得起点可以是任意数量的，目标数组可以是任意数量的。我只处理二维数组，也没有将其扩展到其他维度的意图。我也熟悉pdist、cdist和其他达到我所需解决方案的方法，但我只对einsum感兴趣，因为我想完善我的示例库。任何帮助都将不胜感激。

import numpy as np
origs = np.array([[0.,0.],[1.,0.],[0.,1.],[1.,1.]])
dests = np.asarray([[4.,0.],[1.,1.],[2.,2.],[2.,3.],[0.,5.]])
for i in origs:
    d =np.sqrt(np.einsum("ij,ij->i", i-dests, i-dests))
    print("orig {}...dist: {}".format(i,d))

以下是我正在寻找的结果...

orig [ 0.  0.]...dist: [ 4.          1.41421356  2.82842712  3.60555128  5.        ]
orig [ 1.  0.]...dist: [ 3.          1.          2.23606798  3.16227766  5.09901951]
orig [ 0.  1.]...dist: [ 4.12310563  1.          2.23606798  2.82842712  4.        ]
orig [ 1.  1.]...dist: [ 3.16227766  0.          1.41421356  2.23606798  4.12310563]

- user1121588

值得一提的是，对于未来的读者来说，cdist在性能方面仍然会轻松击败np.einsum（在Divakar的示例中，性能相差一个数量级）。 - ali_m

在我的领域（GIS）中，起点的总数可能不到100个，目的地也不到一千个左右。最终，纯Python / NumPy解决方案消除了安装和维护其他库的需要。虽然它们可能很吸引人且速度快，但有时会混淆幕后实际发生的事情。因此，当数量级仅为几秒甚至十几秒时，这就成为一个非问题。而且由于我是教师，我喜欢让学生看到各种解决方案，而不仅仅是最快的。这有助于促进思维的多样性。 - user1121588

1

很好。关于高效计算欧几里得距离的问题在这里每天都会出现几次，所以我的意图只是指向未来读者可能是最快速的解决方案的方向（特别是因为大多数使用numpy的人也倾向于安装scipy）。 - ali_m

我想所有领域都有它们的差异，NumPy和Python在商业GIS软件中的集成历史非常短，直到最近SciPy、Pandas等才能在软件中使用。人们仍然把后两者视为电影类型和一种动物。 - user1121588

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

如果我理解问题正确，你发布的for循环代码在仅考虑2D数组时看起来很通用。现在，如果你想要一个通用的向量化解决方案，并且只需一次调用np.einsum，你可以将broadcasting引入到其中，就像这样 -

d_all = np.sqrt(np.einsum('ijk->ij',(origs[:,None,:] - dests)**2))

示例运行 -

In [85]: origs = np.array([[0.,0.],[1.,0.],[0.,1.],[1.,1.]])
    ...: dests = np.asarray([[4.,0.],[1.,1.],[2.,2.],[2.,3.],[0.,5.]])
    ...: 

In [86]: for i in origs:
    ...:     d =np.sqrt(np.einsum("ij,ij->i", i-dests, i-dests))
    ...:     print(d)
    ...:     
[ 4.          1.41421356  2.82842712  3.60555128  5.        ]
[ 3.          1.          2.23606798  3.16227766  5.09901951]
[ 4.12310563  1.          2.23606798  2.82842712  4.        ]
[ 3.16227766  0.          1.41421356  2.23606798  4.12310563]

In [87]: np.sqrt(np.einsum('ijk->ij',(origs[:,None,:] - dests)**2))
Out[87]: 
array([[ 4.        ,  1.41421356,  2.82842712,  3.60555128,  5.        ],
       [ 3.        ,  1.        ,  2.23606798,  3.16227766,  5.09901951],
       [ 4.12310563,  1.        ,  2.23606798,  2.82842712,  4.        ],
       [ 3.16227766,  0.        ,  1.41421356,  2.23606798,  4.12310563]])

根据@hpaulj的评论，您也可以使用np.einsum本身执行平方操作，如下所示-

subts = origs[:,None,:] - dests
d_all = np.sqrt(np.einsum('ijk,ijk->ij',subts,subts))

这里有一个运行时测试，用于将其与先前的方法进行比较，该方法在np.einsum之外进行了平方处理 -

In [7]: def all_einsum(origs,dests):
   ...:     subts = origs[:,None,:] - dests
   ...:     return np.sqrt(np.einsum('ijk,ijk->ij',subts,subts))
   ...: 
   ...: def partial_einsum(origs,dests):
   ...:     return np.sqrt(np.einsum('ijk->ij',(origs[:,None,:] - dests)**2))
   ...: 

In [8]: origs = np.random.rand(400,100)

In [9]: dests = np.random.rand(500,100)

In [10]: %timeit all_einsum(origs,dests)
10 loops, best of 3: 139 ms per loop

In [11]: %timeit partial_einsum(origs,dests)
1 loops, best of 3: 251 ms per loop