【问题标题】:Broadcasting with reduction or extension in Numpy在 Numpy 中进行缩减或扩展的广播
【发布时间】:2017-03-26 23:07:03
【问题描述】:

在下面的代码中,我们计算所有给定点对之间的向量大小。为了在 NumPy 中加速这个操作,我们可以使用广播

import numpy as np
points = np.random.rand(10,3)

pair_vectors = points[:,np.newaxis,:] - points[np.newaxis,:,:]
pair_dists = np.linalg.norm(pair_vectors,axis=2).shape

或外部产品迭代

it = np.nditer([points,points,None], flags=['external_loop'], op_axes=[[0,-1,1],[-1,0,1],None])
for a,b,c in it:
    c[...] = b - a
pair_vectors = it.operands[2]
pair_dists = np.linalg.norm(pair_vectors,axis=2)

我的问题是如何使用广播或外部产品迭代来创建一个格式为 10x10x6 的数组,其中最后一个轴包含一对中两个点的坐标(扩展)。并且以一种相关的方式,是否可以直接使用广播或外积迭代来计算对距离,即生成一个 10x10 形式的矩阵,而无需首先计算差异向量(约简)。

为了澄清,以下代码使用慢循环创建所需的矩阵。

pair_coords = np.zeros(10,10,6)
pair_dists = np.zeros(10,10)
for i in range(10):
    for j in range(10):
        pair_coords[i,j,0:3] = points[i,:]
        pair_coords[i,j,3:6] = points[j,:]
        pair_dists[i,j] = np.linalg.norm(points[i,:]-points[j,:])

这是使用外积迭代计算距离(或应用任何其他函数,该函数采用一对中两个点的 6 个坐标并产生一个标量)的失败尝试。

res = np.zeros((10,10))
it = np.nditer([points,points,res], flags=['reduce_ok','external_loop'], op_axes=[[0,-1,1],[-1,0,1],None])
for a,b,c in it: c[...] = np.linalg.norm(b-a)
pair_dists = it.operands[2]

【问题讨论】:

  • pair_vectors 是 10x10x3。我不明白 10x10x6 会包含什么。
  • 这个norm 是:np.sqrt(np.sum(pair_vectors**2,axis=2))。它消除了所有差异;对大小 3 维度求和并取 sqrt。
  • nditer 有点像列表zip。它协调多个数组的迭代。它在处理broadcasting 方面也超越了zip。但是您已经在 c 代码 (cython) 中使用它来获得真正的速度提升。
  • 我最终在 fortran 和 f2py 中使用了循环,它真的很快。我认为 NumPy 是一个适合后处理的环境,但结果证明对简单的循环操作进行矢量化非常复杂。我想 NumPy 在有一种直接的方法来向量化你的函数时很有用。

标签: python numpy array-broadcasting


【解决方案1】:

这是一种以矢量化方式生成这些数组的方法 -

from itertools import product
from scipy.spatial.distance import pdist, squareform

N = points.shape[0]

# Get indices for selecting rows off points array and stacking them 
idx = np.array(list(product(range(N),repeat=2)))
p_coords = np.column_stack((points[idx[:,0]],points[idx[:,1]])).reshape(N,N,6)

# Get the distances for upper triangular elements. 
# Then create a symmetric one for the final dists array.
p_dists = squareform(pdist(points))

this post 中很少讨论其他矢量化方法,所以也可以看看那里!

【讨论】:

  • pdist 在 c 中使用循环,它不允许将周期性边界计入计数。
  • @liberias 抱歉,我们在这里谈论的是什么界限?你能详细说明一下吗?
  • 假设您要计算差异向量,然后对其应用周期性边界,这意味着如果差异向量中的任何坐标超过某个最大框大小(边界),那么您需要减去或添加盒子的大小。这可以通过广播轻松完成,然后使用符号功能和调节。我认为 pdist 不可能,这对我来说毫无用处。但无论如何 pdist 只是用 C 编译的双循环,然后在 Python 中运行,它不能解决使用 NumPy 的问题。
猜你喜欢
  • 2013-07-08
  • 1970-01-01
  • 1970-01-01
  • 2017-12-15
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2013-08-21
相关资源
最近更新 更多