for循环中的纽比点

Question

我在使用的地方有以下代码纽比.dot来加速我的计算。在

u = numpy.zeros((l, l))
wp = numpy.zeros((l,2))
# some code which edits u and wp
for x in range(N):
    wavg = numpy.dot(wp[:, 0], wp[:, 1])
    wp[:, 0] = 1.0/wavg*numpy.dot(u, numpy.multiply(wp[:, 0], wp[:, 1]))

对于小l，最慢的部分是外环。现在我问自己有没有办法摆脱这个循环？在

编辑：用数学术语来说，这个代码应该是这样的

Answer 1

通过预先计算一个项而得到的一个小的改进

u = np.zeros((l, l))
wp = np.zeros((l,2))
# some code which edits u and wp

m = u*wp[:, 1]
for x in range(N):
    wavg = np.dot(wp[:, 0], wp[:, 1])
    wp[:, 0] = 1.0/wavg*np.dot(m, wp[:, 0])

但我们可以做得更好，而不是每次都计算平均值，我们可以在最后一次迭代中完成：

^{pr2}$

但还有一件事我们可以做——这个循环可以用矩阵求幂来代替：

m = u*wp[:, 1]
wp[:, 0] = np.linalg.matrix_power(m, N-1).dot(wp[:, 0])
wavg = np.dot(wp[:, 0], wp[:, 1])
wp[:, 0] = 1.0/wavg*np.dot(m, wp[:, 0])

不幸的是，这似乎要慢得多。但是，如果你能预先计算np.linalg.matrix_power(m, N-1)，那么它会快得多

Answer 2

我仍然对数组的形状以及代码和等式之间的关系感到困惑。在

但只要看看这个方程，我想可以计算为：

In [515]: n,m = 3,4

In [516]: U = np.ones((n,m))

In [517]: w = np.ones((m,))

In [518]: f = np.ones((m,))

In [519]: np.einsum('ij,j,j->i',U,w,f)
Out[519]: array([ 4.,  4.,  4.])

目前，我关心的是尺寸匹配，而不是最终值。这个计算很简单，不需要爱因斯坦符号，但是einsum使得翻译几乎是机械的。在

等效的dot是

^{pr2}$

由于这是随着时间的推移而演变的，迭代的f取决于前一个的值（以及这个外部值w(t)），所以很难消除这个循环；我们可以使内容更快。在