【问题标题】:einsum and distance calculationseinsum 和距离计算
【发布时间】:2015-11-16 05:37:10
【问题描述】:

我已经搜索了一种解决方案,以使用 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]

【问题讨论】:

  • 对于未来的读者来说,值得一提的是,cdist 在性能方面仍将轻松击败np.einsum(在 Divakar 的示例中是一个数量级)。
  • @ali_m 在我的领域 (GIS) 中,起点的总数可能少于 100 个,而目的地的总数可能少于一千个左右。最后,纯 python/numpy 解决方案消除了安装和维护其他库的需要。尽管它们可能很有吸引力且速度很快,但它们有时会混淆幕后实际发生的事情。因此,当数量级仅达到几秒甚至几十秒时,它就变得不成问题了。自从我教书以来,我喜欢让学生看到各种各样的解决方案,而不仅仅是最快的。它促进思维的多样性。
  • 很公平。关于有效计算欧几里得距离的问题往往在这里每天出现好几次,所以我的意图只是将未来的读者指向可能是最快解决方案的方向(特别是因为大多数使用 numpy 的人也倾向于安装 scipy)。
  • 我想所有领域都有其差异,numpy 和 python 被集成到商业 GIS 包中使用的历史很短,而且直到最近,SciPy、Pandas 等才可以在软件。人们仍然将后两者视为电影类型和熊。

标签: python arrays numpy euclidean-distance


【解决方案1】:

如果我正确理解了这个问题,那么在仅考虑 2D 数组时,您发布的 for 循环代码对我来说是通用的。现在,如果您希望通过一次调用 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]])

根据comments by @hpaulj,您也可以使用np.einsum 本身执行平方,就像这样 -

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

这是一个运行时测试,用于将其与之前在 squaring 外部完成 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

【讨论】:

  • einsum 也可以用来取广场。
  • @hpaulj 是通过像这样重复:np.sqrt(np.einsum('ijk,ijk->ij',origs[:,None,:] - dests,(origs[:,None,:] - dests))) 还是有一些更短的代码可以这样做?
  • 是的;虽然不是一行。
  • @hpaulj 不错,它似乎相当快!您想将其发布为答案吗?以我认为的方式分享它会很高兴。
  • 太棒了!我玩过椭圆无济于事,但我没有考虑过没有。我将重新检查文档,但您已经解决了当前的任务。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2012-01-19
  • 2014-08-28
  • 1970-01-01
相关资源
最近更新 更多