加速科学计算的python程序

Question

我有以下Python代码：

def P(z, u0):
    x = np.inner(z, u0)
    tmp = x*u0
    return (z - tmp)


def powerA2(A, u0):
    x0 = np.random.rand(len(A))
    for i in range(ITERATIONS):
        x0 = P(np.dot(A, x0), u0)
        x0 = x0 / np.linalg.norm(x0)
    return (np.inner(np.dot(A, x0), x0))

np 是 numpy 包。

我想在大小为 100,000 * 100,000 的矩阵上运行此代码，但似乎这个程序没有快速运行的机会（我需要运行它很多次，大约 10,000 次）。

有没有什么技巧，如多线程，在这里起作用的机会？

还有其他什么可以帮助加速吗？

- SomeoneHAHA

---

100000的平方在64位实数中将占用大约72GB的内存。这对你来说不是问题吗？我认为加快速度的第一步应该是改用编译语言。我知道numpy正在幕后执行此操作，但我猜测通过在编译语言中本地实现此操作并在Python中调用该例程，可以避免许多内存访问。 - haraldkl

这会是一个很大的改进吗？这个模块是复杂程序的一部分，将其转换为编译语言需要数天时间。 - SomeoneHAHA

看看serge-sans-paille的回答，我不知道pythran，但他展示的基本上是一个3倍的加速。当然这取决于问题，但我猜对于更大的问题，这变得更加重要。正如我所说，你可以考虑用编译语言实现它，并从Python中调用它，虽然这有点麻烦，但应用程序的其余部分可能保持不变。这个Pythran建议基本上为你自动化了这个过程。 - haraldkl

仅使用numpy，我能得到的最佳加速比是1.8。您可以从ITERATIONS循环中删除linalg.norm行以获得适当的加速，但我猜那是为了稳定性而存在的？您可以使用x0 /= math.sqrt(np.dot(x0, x0))而不是x0 = x0/np.linalg.norm(x0)来加快该行的速度。通过使用out=参数并在原地执行一些操作，可以获得更多速度。令人惊讶的是（对我来说），np.einsum没有在任何地方提供帮助； np.inner和np.dot非常快。 - askewchan

Answer 1

你可以考虑使用Pythran。编译以下代码（norm.py）：

#pythran export powerA2(float [][], float[])
import numpy as np

def P(z, u0):
    x = np.inner(z, u0)
    tmp = x*u0
    return (z - tmp)

def norm(x):
    return np.sqrt(np.sum(np.abs(x)**2))

def powerA2(A, u0):
    ITERATIONS = 100
    x0 = np.random.random(len(A))
    for i in range(ITERATIONS):
        x0 = P(np.dot(A, x0), u0)
        x0 = x0 / norm(x0)
    return (np.inner(np.dot(A, x0), x0))

使用：

pythran norm.py

产生以下加速：

$ python -m timeit -s 'import numpy as np; A = np.random.rand(100, 100); B = np.random.random(100); import norm' 'norm.powerA2(A, B)'
100 loops, best of 3: 3.1 msec per loop
$ pythran norm.py -O3 -march=native
$ python -m timeit -s 'import numpy as np; A = np.random.rand(100, 100); B = np.random.random(100); import norm' 'norm.powerA2(A, B)'
1000 loops, best of 3: 937 usec per loop

- serge-sans-paille

---

谢谢！它是否仍然适用于（更）大的矩阵？ - SomeoneHAHA

它将数据存储为简单的平面数组，并尽可能避免使用临时变量，因此只要矩阵适合内存，我会说是的。 - serge-sans-paille

Answer 2

只是想确认一下：您想做10 ^ 4次某个操作，每次操作都是10 ^ 10…即使您的操作是O（1），这仍然是10 ^ 14次操作，这是一个非常困难的问题（正如haraldkl在他的评论中指出的那样，这也会占用大量内存）。只是想确认一下：您要调用powerA2 10,000次，还是10,000是您期望的ITERATIONS值。如果是前者，您可以使用线程（或更好的是，单独的进程）来实现一些并行化，但我不知道这是否足够；如果是后者，除非有我错过的技巧，否则每个循环迭代的输入似乎都不可并行化，因为每个迭代的输入都取决于先前的输出。可能有一种方法可以在GPU上完成这项工作（我希望至少有一种有效的方法可以通过使用矢量化快速地进行大量的规范化处理）。

回应评论后的编辑：cpython（最常见的Python实现）有全局解释器锁（GIL）；一些其他Python实现（jython，ironpython）则没有；参见https://wiki.python.org/moin/GlobalInterpreterLock。

请注意，潜在的阻塞或长时间运行的操作，例如I/O、图像处理和NumPy数值计算，会在GIL之外发生。因此，只有在多线程程序中花费大量时间在GIL内部解释CPython字节码时，GIL才会成为瓶颈。

据我所知，使用numpy进行线程操作应该是可行的，并且不会受到太大的瓶颈，但是除非有一些我不知道的数学技巧，否则你的问题看起来很难转换为线程。

- Foon

---

我读到Python有解释器锁，那么在Python中多个线程可以同时工作吗？ - SomeoneHAHA

Answer 3

通过这种方式重新定义函数，与未编译的serge-sans-paille版本相比，我可以获得10%的改进：

def P0(z, u0):
    x = np.inner(z, u0)
    x *= u0
    return (z - x)

def norm0(x):
    return np.sqrt(np.sum(x*x))

def powerA20(A, u0):
    ITERATIONS = 100
    x0 = np.random.random(len(A))
    for i in range(ITERATIONS):
        x0 = P0(np.dot(A, x0), u0)
        x0 /= norm0(x0)
    return (np.inner(np.dot(A, x0), x0))

使用*= u0代替x = x*u0可以避免在RAM中不必要的变量副本，稍微加快程序速度。 此外，在这种情况下不需要使用abs。最后，x*x比x**2略快。

- Federico Barabas