复数元素的Numpy矩阵向量乘法

Question

矩阵/向量乘法看起来是这样的，工作正常：

A = np.array([[1, 0],
              [0,1]])      # matrix
u = np.array([[3],[9]])    # column vector
U = np.dot(A,u)            # U = A*u
U

如果元素Xx,Xy,Yx,Yy, u,v都是二维或三维数组，有没有办法保持相同的结构？

A = np.array([[Xx, Xy],
              [Yx,Yy]])
u = np.array([[u],[v]])
U = np.dot(A,u)

理想的结果是这样的，而且效果很好（因此我的问题的答案是“很高兴拥有”）：

U = Xx*u + Xy*v
V = Yx*u + Yy*v

成功故事

THX到@loopy walt和@Ivan

1。它可以工作：np.einsum('ijmn,jkmn->ikmn', A, w)和(w.T@A.T).T都再现了正确的逆变坐标变换。

2。次要问题：np.einsum('ijmn,jkmn->ikmn', A, w)和(w.T@A.T).T都返回一个带有一对方括号的数组太多（它应该是二维的）。请参见@Ivan的答案为什么：

array([[[   0,   50,  200],
        [ 450,  800, 1250],
        [1800, 2450, 3200]]])

3。要优化：构建块矩阵非常耗时：

A = np.array([[Xx, Xy],
              [Yx,Yy]])
w = np.array([[u],[v]])

所以总的返回时间是np.einsum('ijmn,jkmn->ikmn', A, w)和(w.T@A.T).T更大的。也许这是可以优化的

- p_einstein np.einsum('ijmn,jkmn->ikmn', A, w)        9.07 µs ± 267 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
- p_tensor: (w.T@A.T).T                                7.89 µs ± 225 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
- p_vect:   [Xx*u1 + Xy*v1, Yx*u1 + Yy*v1]             2.63 µs ± 173 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
- p_vect: including  A = np.array([[Xx, Xy],[Yx,Yy]])  7.66 µs ± 290 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Answer 1

这本质上是一个matrix block运算，可以从简单的矩阵运算中概括出来。您的示例很好，因为您这样定义了块矩阵（我将它们重命名为A_p和u_p，以避免与矩阵u冲突）：

A_p = np.array([[Xx, Xy],
                [Yx, Yy]])
u_p = np.array([[u], [v]])

如果你看一个标准的矩阵运算，就形状而言，你有(r, s)遇到(s, t)，这最终导致了一个形状为(r, t)的矩阵。“s”维度在该过程中减少。如果您还不熟悉Einstein Summation和^{}，那么应该快速查看一下。给定矩阵M和N，矩阵操作如下所示：

out = np.zeros((r, t))
for i in range(r):
   for j in range(s):
      for k in range(t):
         out[i, k] += M[i, j]*N[j, k]

当你意识到操作是如何进行的时，这变得非常直观：覆盖整个多维空间，我们在输出矩阵上累积产品

对于np.einsum，这将非常类似（这里是M=A，这里是N=u）：

>>> np.einsum('ij,jk->ik', A, u)
array([[3],
       [9]])

使用坐标：i、j和k而不是维度大小来思考是很有用的：r、s和t，就像上面使用np.einsum时一样。所以我会在剩下的答案中保留前一个符号

---

现在来看块矩阵，输入矩阵本身包含矩阵块。矩阵M是形状(i, j, (m, n))（即(i, j, m, n)），其中(k, l)表示块的形状。同样，对于N，我们有(j, k, m, n)。我们基本上用2x2块{}替换了标量形状的{}，仅此而已

因此，您可以直接推断A_p*u_p操作为：

>>> np.einsum('ijmn,jkmn->ikmn', A_p, u_p)
array([[[[  0,  50],
         [200, 450]]],


       [[[  0, 110],
         [440, 990]]]])

您可以将其与“手工制作”结果进行比较：

>>> np.stack((Xx*u + Xy*v, Yx*u + Yy*v))[:, None]

注意，这可以用于任何(m, n)块和任何矩阵大小

您现在可以尝试使用类似于A*u的循环来实现A_p*u_p

---

编辑：在@loopy walt的answer上展开，他向您展示了可以通过以下方式获得相同的结果：

>>> (u_p.T @ A_p.T).T

实际上^{}可以处理多维数组，并将以以下方式处理（对于高于2的维度数）：两个输入都将视为驻留在最后两个索引中的矩阵堆栈。输入形状为(i, j, m, n)和(j, k, m, n)。为了使用__matmul__，我们需要将块维度（即(m, n)）拉到第一个位置，并将矩阵维度（分别为(i, j)和(j, k)）放在最后。需要采取以下步骤

• 转置将产生颠倒维度顺序的效果：(n, m, j, i)和(n, m, k, j)

• 然后我们交换两个操作数，我们设法以(*, k, j) x (*, j, i) = (*, k, i)的形式计算一些东西，其中* = (n, m)

• 转换输出会导致(i, k, m, n)的形状

Answer 2

借用@Ivan的A_p和u_p，我想宣传以下更简单的方法，利用__matmul__的内置广播。我们需要做的就是颠倒尺寸的顺序：

(u_p.T@A_p.T).T

这给出了与Ivan的einsum相同的结果，但更灵活一些，可以提供完整的广播，而无需每次都计算出相应的格式字符串

np.allclose(np.einsum('ijmn,jkmn->ikmn', A_p, u_p),(u_p.T@A_p.T).T)
# True