我有一个k*n 矩阵X 和一个k*k 矩阵A。对于X 的每一列，我想计算标量
X[:, i].T.dot(A).dot(X[:, i])
（或者，在数学上，Xi' * A * Xi）。
目前，我有一个for 循环：
out = np.empty((n,))
for i in xrange(n):
out[i] = X[:, i].T.dot(A).dot(X[:, i])
但由于 n 很大，如果可能的话，我想更快地执行此操作（即使用一些 NumPy 函数而不是循环）。
【问题讨论】：
标签： numpy
这似乎做得很好：
(X.T.dot(A)*X.T).sum(axis=1)
编辑：这有点快。 np.einsum('...i,...i->...', X.T.dot(A), X.T)。如果 X 和 A 是 Fortran 连续的，则两者都可以更好地工作。
【讨论】：
n=10000, k=10，我的代码是76.2ms，新代码是1.64ms。不错！
您可以使用numpy.einsum:
np.einsum('ji,jk,ki->i',x,a,x)
这将得到相同的结果。让我们看看它是否更快：
看起来dot 仍然是最快的选择，特别是因为它使用线程 BLAS，而不是在一个内核上运行的einsum。
import perfplot
import numpy as np
def setup(n):
k = n
x = np.random.random((k, n))
A = np.random.random((k, k))
return x, A
def loop(data):
x, A = data
n = x.shape[1]
out = np.empty(n)
for i in range(n):
out[i] = x[:, i].T.dot(A).dot(x[:, i])
return out
def einsum(data):
x, A = data
return np.einsum('ji,jk,ki->i', x, A, x)
def dot(data):
x, A = data
return (x.T.dot(A)*x.T).sum(axis=1)
perfplot.show(
setup=setup,
kernels=[loop, einsum, dot],
n_range=[2**k for k in range(10)],
logx=True,
logy=True,
xlabel='n, k'
)
【讨论】：
for 循环更快...值得一试
for 循环 cython/numpy for 循环无关紧要。时间真的不在循环中。
n=10000，这优于我的原始代码（76.2ms vs. 1.48ms）。
einsum 链接；很高兴知道它可以做什么。太糟糕了，这不是最快的解决方案。 （非常令人惊讶的是，它比 @IanH 在一个内核上的 n=10000, k=10 解决方案更快）
除非您将整个事情并行化，否则您无法更快地做到这一点：每列一个线程。你仍然会使用循环——你无法摆脱它。
Map reduce 是解决这个问题的好方法：map step multiples，reduce step sums。
【讨论】：