Python 3：优化 scipy 数组的求和答案

【问题标题】：Python 3: Optimizing summation over scipy arraysPython 3：优化 scipy 数组的求和
【发布时间】：2013-09-26 07:18:18
【问题描述】：

我目前正在解决一个问题，我必须对 scipy/numpy 数组的特定条目进行求和，并且我正在寻找一种方法来完全摆脱所有 Python for 循环。我在 Mac OS X 上使用 Python 3.3。下面是我正在做的总和之一的示例，为了示例的缘故，我对填充有随机整数的数组的条目进行求和。

from scipy import ones, conjugate, sum, random

n = 5
M = random.randint(5,size=(4*n**2,4*n**2))
H = sum((M[i+1,:2*n**2]*M[i,:2*n**2].conjugate()).sum() * (-M[i,:2*n**2]*M[i+1,:2*n**2].conjugate()).sum() for i in range(0,2*n**2,2))

我首先计算两个矩阵条目的乘积，然后对一半的列求和。我这样做了两次，然后以两步对一半以上的行求和。

这可能看起来很奇怪，但我正在使用哈密顿量来构建格子上的系统，其中每一行对应于某个格点，偶数行和奇数行代表向上或向下旋转。 n 最终会很大，我需要加快总结速度。

现在，我无法弄清楚如何摆脱 for 循环。我尝试在行索引中使用 range() 参数来执行此操作，但这并没有得到相同的结果。

谢谢！

【问题讨论】：

只有我一个人，还是你只使用了数组的四分之一？
是的，我通过对角化矩阵得到数组，M 对应于填充了特征向量的矩阵。我已经转换了总和，使得循环最短。其他总和将使用数组的不同部分。

标签： python arrays optimization numpy scipy

【解决方案1】：

在我看来你可以做到

H_bis = np.sum(M[1:2*n**2:2, :2*n**2] * M[:2*n**2:2, :2*n**2].conjugate(), axis=1)
H_bis = H_bis * H_bis.conjugate()
H_bis = -np.sum(H_bis)

【讨论】：

【解决方案2】：

这个怎么样：

# We need a complex test case to make sure the conjugate works properly
M = (random.randint(5,size=(4*n**2,4*n**2)) +
     random.randint(5,size=(4*n**2,4*n**2))*1j)

H_bis = np.sum(M[1:2*n**2:2, :2*n**2] * M[:2*n**2:2, :2*n**2].conjugate(),
               axis=-1)
H_bis *= np.sum(-M[:2*n**2:2, :2*n**2] * M[1:2*n**2:2, :2*n**2].conjugate(),
                axis=-1)
H_bis = np.sum(H_bis)
>>> np.allclose(H, H_bis)
True

【讨论】：

我想你可能错过了那里的共轭。
你说得对，我看错了索引，它之所以有效，是因为他的数据是真实的......
只是我还是这个比原来慢？
对于n = 5 的 OP 值，我得到了 x8.5 的加速。但是对于更大的n 事情甚至可以解决，并且在我的系统上，对于n = 20，OP 的方法实际上要快一点。可能是因为这种矢量化方法，尤其是在没有优化的情况下，会占用更多内存。