<pre><code>X = np.einsum('ik,jk->ij', np.conj(x), x)
</code></pre>
<p>相当于</p>
^{pr2}$
<p><a href="http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.einsum.html" rel="nofollow">^{<cd1>}</a>
取产品总数。下标<code>'ik,jk->ij'</code>告诉<code>np.einsum</code>第二个参数,
<code>np.conj(x)</code>是一个下标为<code>ik</code>的数组,第三个参数<code>x</code>有
下标<code>jk</code>。因此,计算所有的乘积<code>np.conj(x)[i,k]*x[j,k]</code>
<code>i</code>,<code>j</code>,<code>k</code>。求和取重复下标<code>k</code>,因此
保留<code>i</code>和<code>j</code>,它们将成为结果数组的下标。在</p>
<hr/>
<p>例如</p>
<pre><code>import numpy as np
N, M = 10, 20
a = np.random.random((N,M))
b = np.random.random((N,M))
x = a + b*1j
def orig(x):
X = np.empty((x.shape[0], x.shape[0]), dtype='complex128')
for i in range(x.shape[0]):
for j in range(x.shape[0]):
X[i, j] = np.vdot(x[i], x[j])
return X
def alt(x):
return np.einsum('ik,jk->ij', np.conj(x), x)
assert np.allclose(orig(x), alt(x))
</code></pre>
<hr/>
<pre><code>In [307]: %timeit orig(x)
10000 loops, best of 3: 143 µs per loop
In [308]: %timeit alt(x)
100000 loops, best of 3: 8.63 µs per loop
</code></pre>