使用 numpy 在没有 for 循环的情况下更新节点值答案

【问题标题】：Update nodal values without for loop using numpy使用 numpy 在没有 for 循环的情况下更新节点值
【发布时间】：2019-04-01 16:09:35
【问题描述】：

我正在尝试根据元素值更新 mesh 上的节点值。

在一个数组faces 中，我定义了一个元素节点的ID（假设我只有两个元素）：

faces = np.array([[0, 1, 2], [1, 3, 2]])

数组force_el 包含，比方说，作用在元素每个节点上的力：

force_el = np.array([[0.7, 1.1], [1.2, 0.3]])

现在我想更新节点力force_node：

force_node = np.zeros((4, force_el.shape[1]))
for face, fel in zip(faces, force_el):
    force_node[face.ravel(), :] += fel

所以结果是：

>>> force_node
array([[0.7, 1.1],
       [1.9, 1.4],
       [1.9, 1.4],
       [1.2, 0.3]])

由于此更新必须进行多次（大约 100k-1m 次），我正在尝试对其进行优化，但找不到好的解决方案。

【问题讨论】：

标签： python performance numpy vectorization

【解决方案1】：

你可以使用一些matrix-multiplicationforce -

out_nrows = 4 # number of nodes
mask = np.zeros((len(faces),out_nrows),dtype=bool)
np.put_along_axis(mask,faces,True,axis=1)
force_node_out = mask.T.dot(force_el)

force_el 中的列数较少，我们还可以使用 np.bincount 以获得更好的性能 -

out_nrows = 4 # number of nodes
out = np.zeros((out_nrows, force_el.shape[1]))
n = faces.shape[1]
l = force_el.shape[1]
for i in range(n):
    for j in range(l):
        out[:,j] += np.bincount(faces[:,i],force_el[:,j],minlength=out_nrows)

时间安排 -

In [35]: # Setup data (from OP's comments)
    ...: np.random.seed(0)
    ...: faces=np.array([np.random.choice(1800,3,replace=0) for i in range(3500)])
    ...: force_el = np.random.rand(len(faces),3)

In [36]: %%timeit # Original loopy soln
    ...: out_nrows = 1800
    ...: force_node = np.zeros((out_nrows, force_el.shape[1]))
    ...: for face, fel in zip(faces, force_el):
    ...:     force_node[face.ravel(), :] += fel
100 loops, best of 3: 16.1 ms per loop

In [37]: %%timeit # @RafaelC's soln with np.add.at
    ...: force_node = np.zeros((1800, force_el.shape[1]))
    ...: np.add.at(force_node, faces, force_el[:,None])
100 loops, best of 3: 2.45 ms per loop

In [38]: %%timeit # Posted in this post that uses matrix-multiplication
    ...: out_nrows = 1800
    ...: mask = np.zeros((len(faces),out_nrows),dtype=bool)
    ...: np.put_along_axis(mask,faces,True,axis=1)
    ...: force_node_out = mask.T.dot(force_el)
10 loops, best of 3: 38.4 ms per loop

In [39]: %%timeit # Posted in this post that uses bincount
    ...: out_nrows = 1800
    ...: out = np.zeros((out_nrows, force_el.shape[1]))
    ...: n = faces.shape[1]
    ...: l = force_el.shape[1]
    ...: for i in range(n):
    ...:     for j in range(l):
    ...:         out[:,j]+=np.bincount(faces[:,i],force_el[:,j],minlength=out_nrows)
10000 loops, best of 3: 149 µs per loop

【讨论】：

感谢分享您的解决方案。在我的帖子中，我没有提到系统远大于 4 个节点。通过测试您的解决方案和@RafaelC 假设 1800 个节点和 10000 个时间步长的解决方案，我得到了以下结果： - Loopy 解决方案：~159 s - 您的解决方案：~231s - RafaelC 的解决方案：~23 s 所以我猜对于大型系统 rafael 的方式是要带的那个。
@David 使用您的实际数据集 - force_el 的脸的形状是什么？
两个数组的形状都是 (3500, 3)，而 force_node.shape == (1800, 3)。抱歉，系统的大小不清楚。只是为了给你更多的背景知识，我正在用三角形元素建模一个空心圆柱体。因此，面数与节点数相关，反之亦然。
@David 为此类数据集添加了计时。看一看？你最后有没有给bincoount一个时间？
这是一个巨大的性能提升，谢谢！使用与上述完全相同的配置，您的 bincount 解决方案只需约 1.5 秒。非常感谢。

【解决方案2】：

您应该尽可能利用numpy 广播。

使用np.add.at

np.add.at(force_node, faces, force_el[:,None])

【讨论】：

谢谢@RafaelC！这实际上是我正在寻找的方法。您的答案可能已被标记为已接受的答案，我相信这将对我/将来对这个问题的其他用户有所帮助。由于我真的很关心性能，@Divakar 的解决方案是我的具体情况的最佳解决方案。
@David 无论如何！很高兴我能帮上忙。 np.add.at 的优点是可读性强且直观，但性能方面 Divakar 的矩阵 mult 要好得多。无论如何，迪瓦卡是神，应该接受他的回答哈哈。良好的编码：}