【发布时间】：2022-01-05 07:39:18
【问题描述】：
我正在尝试确定 Python/Numpy 是否是开发我的数字软件的可行替代方案，该软件已经在 C++ 中可用。为了在 Python/Numpy 中获得性能，需要“矢量化”代码。但事实证明，一旦我离开非常简单的示例，我就很难对代码进行矢量化（我不是在谈论 SIMD 指令，而是在谈论没有循环的“高效 Numpy 代码”）。这是我想在 Python/Numpy 中高效实现的算法。
- 创建一个 numpy 数组，其中包含：1.0、1.0 + 1/n、1.0 + 2/n、...、2.0
- 对于数组中的每个 u，使用牛顿法计算 x^2 - u 的根，当 |dx| 时停止
- 对结果数组的所有元素求和
这是我要加速的 Python 算法
import numpy as np
n = 1000000
data = np.arange(1.0, 2.0, 1.0 / n)
def newton(u):
x = 2.0
while True:
f = x**2 - u
df_dx = 2 * x
dx = f / df_dx
if (abs(dx) <= 1.0e-7):
break
x -= dx
return x
result = map(newton, data)
print result[n - 1]
这里是 C++11 中算法的一个版本
#include <iostream>
#include <vector>
#include <cmath>
int main (int argc, char const *argv[]) {
auto n = std::size_t{100000000};
auto v = std::vector<double>(n + 1);
for(size_t k = 0; k < v.size(); ++k) {
v[k] = 1.0 + static_cast<double>(k) / n;
}
auto result = std::vector<double>(n + 1);
for(size_t k = 0; k < v.size(); ++k) {
auto x = double{2.0};
while(true) {
auto f = double{x * x - v[k]};
auto df_dx = double{2 * x};
auto dx = double{f / df_dx};
if (std::abs(dx) <= 1.0e-7) {
break;
}
x -= dx;
}
result[k] = x;
}
auto somme = double{0.0};
for(size_t k = 0; k < result.size(); ++k) {
somme += result[k];
}
std::cout << somme << std::endl;
return 0;
}
在我的机器上运行需要 2.9 秒。有没有办法制作一个快速的 Python/Numpy 算法来做同样的事情（我愿意得到慢不到 5 倍的东西）。
谢谢。
【问题讨论】：
-
numpy 的好处是它在幕后进行矢量化，消除了很多（但不是全部）你在编写 C++ 时需要做的优化。当然，您肯定需要使用 numpy 的所有工具来获得最高性能的代码，例如使用
ndarrays 代替列表列表。 -
@MattDMo：我正在寻找一个实际有效的代码。我只是不知道该怎么做，这就是我寻求帮助的原因。
-
'向量化'对于串行操作很困难 - 步骤
i中的计算取决于步骤i-1的结果。当所有步骤都可以并行完成时会更容易 - 即评估顺序无关紧要。 -
@hpaulj：这里不是这样，因为 newton(u1) 和 newton(u2) 是完全独立的。只是在 Python/Numpy 中没有办法有效地做到这一点。我发现的唯一解决方案是使用 Cython 或 Numba，这对我来说是不行的。我不想为这么简单的事情“破解”。
-
在
newton循环中，x一步依赖于上一步的x（x -= dx和dx本身是x的函数）。这就是我所说的串行或迭代解决方案的意思。