numpy - einsum 表示法：矩阵堆栈与向量堆栈的点积

Question

我想将 m*m 矩阵的 n-dim 堆栈乘以向量的 n-dim 堆栈（长度为 m），以便生成的 m*n 数组包含矩阵和向量的点积的结果在第 n 个条目中：

vec1=np.array([0,0.5,1,0.5]); vec2=np.array([2,0.5,1,0.5])
vec=np.transpose(n.stack((vec1,vec2)))
mat = np.moveaxis(n.array([[[0,1,2,3],[0,1,2,3],[0,1,2,3],[0,1,2,3]],[[-1,2.,0,1.],[0,0,-1,2.],[0,1,-1,2.],[1,0.1,1,1]]]),0,2)
outvec=np.zeros((4,2))
for i in range(2):
    outvec[:,i]=np.dot(mat[:,:,i],vec[:,i])

受Element wise dot product of matrices and vectors这篇文章的启发，我尝试了 einsum 中索引组合的所有不同扰动，并发现

np.einsum('ijk,jk->ik',mat,vec)

给出正确的结果。

不幸的是，我真的不明白这一点 - 我假设我在 'ijk,jk' 部分重复输入 k 的事实意味着我将 k 相乘并求和。我已经尝试阅读文档https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.einsum.html，但我仍然不明白。

（包括我之前的尝试，

 np.einsum('ijk,il->ik', mat, vec)

我什至不确定这意味着什么。当我删除它时，索引 l 会发生什么？）

提前致谢！

Answer 1

阅读Einstein summation notation。

基本上，规则是：

没有->

• 输入中重复的任何字母都表示要进行乘法和求和的轴
• 输入中未重复的任何字母都包含在输出中

带有->

• 输入中重复的任何字母都表示要被倍增的轴
• 任何不在输出中的字母都代表要求和的轴

因此，例如，矩阵 A 和 B 具有相同的形状：

np.einsum('ij, ij',       A, B)  # is A ddot B,                returns 0d scalar
np.einsum('ij, jk',       A, B)  # is A dot  B,                returns 2d tensor
np.einsum('ij, kl',       A, B)  # is outer(A, B),             returns 4d tensor
np.einsum('ji, jk, kl',   A, B)  # is A.T @ B @ A,             returns 2d tensor
np.einsum('ij, ij -> ij', A, B)  # is A * B,                   returns 2d tensor
np.einsum('ij, ij -> i' , A, A)  # is norm(A, axis = 1),       returns 1d tensor
np.einsum('ii'             , A)  # is tr(A),                   returns 0d scalar

Answer 2

In [321]: vec1=np.array([0,0.5,1,0.5]); vec2=np.array([2,0.5,1,0.5])
     ...: vec=np.transpose(np.stack((vec1,vec2)))
In [322]: vec1.shape
Out[322]: (4,)
In [323]: vec.shape
Out[323]: (4, 2)

stack 函数的一个好处是我们可以指定一个轴，跳过转置：

In [324]: np.stack((vec1,vec2), axis=1).shape
Out[324]: (4, 2)

为什么要混合使用 np. 和 n.？ NameError: name 'n' is not defined。这种事情差点把我吓跑了。

In [326]: mat = np.moveaxis(np.array([[[0,1,2,3],[0,1,2,3],[0,1,2,3],[0,1,2,3]],[[-1,2.,0
     ...: ,1.],[0,0,-1,2.],[0,1,-1,2.],[1,0.1,1,1]]]),0,2)
In [327]: mat.shape
Out[327]: (4, 4, 2)

In [328]: outvec=np.zeros((4,2))
     ...: for i in range(2):
     ...:     outvec[:,i]=np.dot(mat[:,:,i],vec[:,i])
     ...:     
In [329]: outvec
Out[329]: 
array([[ 4.  , -0.5 ],
       [ 4.  ,  0.  ],
       [ 4.  ,  0.5 ],
       [ 4.  ,  3.55]])

In [330]: # (4,4,2) (4,2)   'kji,ji->ki'

从您的循环中，i 轴（大小 2）的位置很清楚 - 在所有 3 个数组中的最后一个。这为vec 留下了一个轴，我们称之为j。它与最后一个配对（在mat 的i 旁边）。 k 从 mat 延续到 outvec。

In [331]: np.einsum('kji,ji->ki', mat, vec)
Out[331]: 
array([[ 4.  , -0.5 ],
       [ 4.  ,  0.  ],
       [ 4.  ,  0.5 ],
       [ 4.  ,  3.55]])

einsum 字符串通常会自行写入。例如，如果mat 被描述为 (m,n,k)，vec 被描述为 (n,k)，结果为 (m,k)

在这种情况下，只有 j 维度被求和 - 它出现在左侧，但在右侧。最后一个维度，i 在我的符号中，没有求和，因为 if 出现在两边，就像它在你的迭代中一样。我认为这是“随心所欲”。它不是dot 产品的积极组成部分。

实际上，你是在最后一个维度上堆叠，尺寸为 2。通常我们堆叠在第一个，但你将两者都调换到最后。

---

您的“失败”尝试运行，并且可以重现为：

In [332]: np.einsum('ijk,il->ik', mat, vec)
Out[332]: 
array([[12. ,  4. ],
       [ 6. ,  1. ],
       [12. ,  4. ],
       [ 6. ,  3.1]])
In [333]: mat.sum(axis=1)*vec.sum(axis=1)[:,None]
Out[333]: 
array([[12. ,  4. ],
       [ 6. ,  1. ],
       [12. ,  4. ],
       [ 6. ,  3.1]])

j 和 l 维度未出现在右侧，因此将它们相加。它们可以在相乘之前求和，因为它们每个只出现在一个术语中。我添加了None 以启用广播（将ik 与i 相乘）。

np.einsum('ik,i->ik', mat.sum(axis=1), vec.sum(axis=1))

---

如果您堆叠在第一个上，并为 vec (2,4,1) 添加了一个维度，那么它会在 matmul 上添加一个 (2,4,4) 垫子。 mat @ vec[...,None].

In [337]: m1 = mat.transpose(2,0,1)
In [338]: m1@v1[...,None]
Out[338]: 
array([[[ 4.  ],
        [ 4.  ],
        [ 4.  ],
        [ 4.  ]],

       [[-0.5 ],
        [ 0.  ],
        [ 0.5 ],
        [ 3.55]]])
In [339]: _.shape
Out[339]: (2, 4, 1)

Answer 3

einsum 很简单（当您玩过索引排列一段时间后，就是...）。

让我们使用一些简单的东西，一个 2×2 矩阵的三重堆栈和 2×, 数组的三重堆栈

import numpy as np

a = np.arange(3*2*2).reshape((3,2,2))
b = np.arange(3*2).reshape((3,2))

我们需要知道我们将使用einsum 计算什么

In [101]: for i in range(3): 
     ...:     print(a[i]@b[i])                                                                            
[1 3]
[23 33]
[77 95]

我们做了什么？我们有一个索引i，当我们在一个堆叠矩阵和一个堆叠向量（均由i索引）之间执行点积时，该索引是固定的，并且单独的输出行意味着对最后一个索引的求和堆叠矩阵和堆叠向量的唯一索引。

这很容易在 einsum 指令中编码

• 我们希望使用相同的i 索引来指定矩阵、向量以及输出，
• 我们想沿着最后一个矩阵索引和剩余的向量索引减少，比如k
• 我们希望输出中的列数与每个堆叠矩阵中的行数一样多，例如j

因此

In [102]: np.einsum('ijk,ik->ij', a, b)                                                                   
Out[102]: 
array([[ 1,  3],
       [23, 33],
       [77, 95]])

我希望我对如何获得正确指令的讨论清晰、正确且有用。