计算 - numpy python bug答案

【问题标题】：Calculation - numpy python bug计算 - numpy python bug
【发布时间】：2011-07-20 16:09:23
【问题描述】：

我正在使用 NumPy 通过大盒子和小盒子之间的 Aperture 来计算 Y 截距。我在大盒子里有超过 100.000 个粒子，在小盒子里有大约 1000 个。这样做需要很多时间。所有 self.YD、self.XD 都是我要相乘的非常大的数组。

PS：ind 是需要相乘的值的索引。在我的代码中该行之前我有一个非零条件。

有什么想法可以用更简单的方式进行计算吗？

YD_zero = self.oldYD[ind] - ((self.oldYD[ind]-self.YD[ind]) * self.oldXD[ind])/(self.oldXD[ind]-self.XD[ind])

谢谢！

更新

将使用 Numpy 的乘法、除法、减法和所有这些东西。让它更快？或者，如果我拆分计算。例如。

首先这样做：

    YD_zero = self.oldYD[ind] - ((self.oldYD[ind]-self.YD[ind])*self.oldXD[ind])

然后下一行将是：

    YD_zero /= (self.oldXD[ind]-self.XD[ind])

有什么建议吗？！

更新 2

一段时间以来，我一直在尝试解决这个问题，但进展不大。我担心的是分母：

    self.oldXL[ind]-self.XL[ind] == 0

我得到了一些奇怪的结果。

另一件事是非零函数。我已经测试了一段时间了。谁能告诉我这和在 Matlab 中的 find 几乎一样

【问题讨论】：

如果您尝试修剪点集，精度损失是否显着？
我需要非常准确的数据。
更简单的意思是CPU周期更少？
我在 while 循环中有这个，并且一次运行至少 16000 次。所以，无论怎样能让这个计算更快。

标签： python matlab numpy bug-tracking

【解决方案1】：

也许我弄错了，但在 Numpy 中，您可以执行矢量化计算。删除封闭的while 循环并运行它...

YD_zero = self.oldYD - ((self.oldYD - self.YD) * self.oldXD) / (self.oldXD - self.XD)

应该会快很多。

更新：使用 Newton-Raphson 方法迭代求根 ...

unconverged_mask = np.abs(f(y_vals)) > CONVERGENCE_VALUE:
while np.any(unconverged_mask):
    y_vals[unconverged_mask] = y_vals[unconverged_mask] - f(y_vals[unconverged_mask]) / f_prime(y_vals[unconverged_mask])
    unconverged_mask = np.abs(f(y_vals)) > CONVERGENCE_VALUE:

此代码只是说明性的，但它显示了如何使用矢量化代码将迭代过程应用于任何函数f，您可以找到f_prime 的导数。 unconverged_mask 表示当前迭代的结果将只应用于那些尚未收敛的值。

请注意，在这种情况下不需要迭代，Newton-Raphson 将在第一次迭代中给你正确的答案，因为我们正在处理直线。你所拥有的是一个精确的解决方案。

第二次更新

好的，所以您没有使用 Newton-Raphson。要一次性计算YD_zero（y 截距），您可以使用，

YD_zero = YD + (XD - X0) * dYdX

其中dYdX 是渐变，在您的情况下似乎是，

dYdX = (YD - oldYD) / (XD - oldXD)

我假设 XD 和 YD 是粒子的当前 x,y 值，oldXD 和 oldYD 是粒子之前的 x,y 值，X0 是粒子的 x 值光圈。

仍然不完全清楚为什么必须遍历所有粒子，Numpy 可以一次对所有粒子进行计算。

【讨论】：

如果你对它进行矢量化处理，如果你真的需要速度，它还允许稍后移动到 GPU。
我确实对它进行了矢量化，但我还有很多其他的东西在循环中。我认为这比其他任何事情都需要更长的时间。因为我在没有这个的情况下运行它，而且速度要快得多。
我相信循环不在索引上，而是针对迭代牛顿算法。
啊，对了，所以您正在使用 Newton-Raphson 方法在函数中查找零点？如果是这样，更新的代码可能会有所帮助。
另外请注意，由于我们处理的是直线，因此无需多次应用 Newton-Raphson 方程，它会在第一次迭代时正确。

【解决方案2】：

由于所有计算都是按元素完成的，因此在Cython 中重写表达式应该很容易。这将避免在您执行oldYD-YD 等操作时创建的所有那些非常大的临时数组。

另一种可能是numexpr。

【讨论】：

如何将两个数组与 Numexpr 或 Cython 相乘？！您更愿意推荐 Cython 还是 Numexpr？
@theSun：乘法的具体原因是什么？至于Cython vs numexpr，后者应该更容易做对。但是，我个人只使用了前者，所以numexpr 周围可能存在一些我不知道的细微之处。

【解决方案3】：

我肯定会选择numexpr。我不确定numexpr 是否可以处理索引，但我敢打赌以下（或类似的）会起作用：

import numexpr as ne

yold = self.oldYD[ind]
y = self.YD[ind]
xold = self.oldXD[ind]
x = self.XD[ind]
YD_zero = ne.evaluate("yold - ((yold - y) * xold)/(xold - x)")

【讨论】：