加速科学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包。

我对运行100000*100000矩阵的代码很感兴趣，但是这个程序似乎没有机会运行得快（我需要运行它很多次，大约10000次）。

这里有没有可能使用多线程之类的技巧？

还有什么能帮助加速的吗？

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))

有：

^{pr2}$

产生以下加速：

$ 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

Answer 2

只是检查一下：你想对10^10做10^4个操作。。。因此，即使你的操作是O（1），那仍然是10^14个操作，这是一个相当困难的问题（正如haraldkl在评论中指出的，这也在消耗大量内存）只是为了检查：你是要调用powerA2 10000次，还是10000次是你期望的迭代值。如果是前者，您可以使用线程（或者更好的是，单独的进程）来获得一些并行化，但我不知道这是否足够；如果后者，除非我遗漏了一个技巧，否则您的输入看起来不是可并行的，因为每个循环迭代的输入依赖于前一个循环的输出。在GPU上可能有一种方法可以做到这一点（我想至少有一种有效的方法来完成规范化位，这样它就可以通过使用向量化快速地完成大量的工作）

编辑以响应注释：cpython（这是最常见的python实现）有一个全局Interpeter锁（GIL）；其他一些python实现（jython、ironpython）没有；per https://wiki.python.org/moin/GlobalInterpreterLock。在

Note that potentially blocking or long-running operations, such as I/O, image processing, and NumPy number crunching, happen outside the GIL. Therefore it is only in multithreaded programs that spend a lot of time inside the GIL, interpreting CPython bytecode, that the GIL becomes a bottleneck.

据我所知，它应该可以使用线程与numpy和不可怕的瓶颈，但你的问题看起来仍然很难转换成线程，除非有一些数学我遗漏了。在

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而不是{}这样的操作可以避免RAM中不必要的变量副本，从而使程序稍微加快一点。 你也不需要这样。最后，x*x比x**2稍快。在