二维数组每列的外积形成三维数组 - NumPy

Question

假设X是一个M x N的矩阵，用xi表示X的第i列。我想创建一个由M x M矩阵构成的三维N x M x M数组，每个矩阵是xi.dot(xi.T)。

如何使用numpy最优雅地实现？是否可能只使用矩阵运算而不需要循环来完成这个操作？

- Virtuoz

---

听起来像是 np.einsum 的工作。 - user2357112

似乎在这里混淆了（行，列）...使用具体数字，假设我们从（M行，N列）=（2,3）开始。 OP正在使用列向量操作 Ni * Ni.T。每列的长度为Mrows = 2，因此乘积产生一个（2，2）矩阵，即（Mrows，Mrows）。对于每一列进行这样的操作，并将结果堆叠在第三维中，应该得到一个（M x M x N）矩阵。 - Matt P

@MattP OP 可能是沿着第一个轴进行堆叠/存储：for i in range(N): out[i] = X[:,i,None].dot(X[None,:,i])。 - Divakar

@Divakar，这里的广播/索引操作做得不错！今天我在这里学到了新东西。我的解决方案（没有那么花哨，循环可能更多一些，而且可能性能更低），可以得到相同的结果:Result = np.dstack( X[:,i].reshape((nrows,1)) * X[:,i] for i in range(ncols) ). - Matt P

@MattP 更简单的版本是：xi = X[:,i,None]，然后使用OP的代码：xi.dot(xi.T)。 - Divakar

Answer 1

采用 broadcasting 的一种方式 -

X.T[:,:,None]*X.T[:,None]

另一个使用 broadcasting 并在其后交换轴的方法 -

(X[:,None,:]*X).swapaxes(0,2)

另一个涉及广播和多维转置的操作 -

(X[:,None,:]*X).T

使用np.einsum的另一个方法，如果您从循环代码进行转换，则可能更加直观地考虑涉及的迭代器-

np.einsum('ij,kj->jik',X,X)

所有这些方法的基本思想都是我们将最后一个轴展开，使其与其他轴对齐进行逐元素相乘。我们通过将 X 扩展为两个 3D 数组版本来实现这个“互相对抗”的过程。

- Divakar