理解 NumPy 的 einsum
- 2024-12-11 08:47:00
- admin 原创
- 142
问题描述：
怎么np.einsum運作？
给定数组A和B，它们的矩阵乘法然后是转置使用来计算(A @ B).T，或者等效地使用：
np.einsum("ij, jk -> ki", A, B)
解决方案 1：
（注意：这个答案基于我之前写的一篇简短的博客文章。）einsum
做什么einsum？
假设我们有两个多维数组，A和B。现在假设我们想要...
A以特定的方式繁殖B，创造出一系列新的产品；然后也许沿特定轴对这个新数组求和；然后也许
按照特定顺序转置新数组的轴。
这很有可能einsum帮助我们比 NumPy 函数组合（如multiply、sum和）更快、更节省内存地完成此操作transpose。
怎么einsum運作？
这是一个简单（但并非完全无足轻重）的例子。以以下两个数组为例：
A = np.array([0, 1, 2])
B = np.array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
我们将逐个元素相乘A，B然后沿着新数组的行求和。在“普通” NumPy 中，我们会这样写：
>>> (A[:, np.newaxis] * B).sum(axis=1)
array([ 0, 22, 76])
因此，这里的索引操作将A两个数组的第一个轴对齐，以便可以广播乘法。然后将乘积数组的行相加以返回答案。
现在如果我们想使用einsum，我们可以写：
>>> np.einsum('i,ij->i', A, B)
array([ 0, 22, 76])
签名字符串'i,ij->i'是这里的关键，需要稍微解释一下。你可以把它想象成两半。在左侧（ 的左侧） ，->我们标记了两个输入数组。在 的右侧->，我们标记了我们想要的最终数组。
接下来会发生什么：
A有一个轴；我们将其标记为i。并且B有两个轴；我们将轴 0 标记为i，将轴 1 标记为j。通过在两个输入数组中重复标签
i，我们告诉einsum这两个轴应该相乘。换句话说，我们将数组A与数组的每一列相乘B，就像这样A[:, np.newaxis] * B做一样。请注意，
j在我们期望的输出中，不会出现标签；我们刚刚使用了i（我们希望最终得到一个一维数组）。通过省略标签，我们告诉einsum沿此轴求和。换句话说，我们正在对产品的行求和，就像这样.sum(axis=1)做一样。
这基本上就是您需要了解的有关使用 的全部内容einsum。稍微玩一下会有所帮助；如果我们将两个标签都留在输出中，'i,ij->ij'，我们将返回一个 2D 产品数组（与 相同A[:, np.newaxis] * B）。如果我们说没有输出标签，'i,ij->，我们将返回一个数字（与 相同(A[:, np.newaxis] * B).sum()）。
然而，最棒的einsum是，它不会先构建一个临时的乘积数组；它只是在计算过程中对乘积进行求和。这可以大大节省内存使用量。
一个稍微大一点的例子
为了解释点积，这里有两个新数组：
A = array([[1, 1, 1],
[2, 2, 2],
[5, 5, 5]])
B = array([[0, 1, 0],
[1, 1, 0],
[1, 1, 1]])
我们将使用 计算点积。下图显示了和np.einsum('ij,jk->ik', A, B)的标签以及我们从函数中获得的输出数组：A`B`
您可以看到标签j重复了 - 这意味着我们将的行A与的列相乘B。此外，标签j不包含在输出中 - 我们正在将这些乘积相加。标签i和k保留用于输出，因此我们返回一个二维数组。
j将此结果与未对标签求和的数组进行比较可能会更清楚。下面左侧您可以看到写入后生成的 3D 数组np.einsum('ij,jk->ijk', A, B)（即我们保留了标签j）：
求和轴j给出预期的点积，如右图所示。
一些练习
为了更好地理解einsum，使用下标符号实现熟悉的 NumPy 数组操作会很有用。任何涉及乘法和求和轴组合的内容都可以使用 来编写 einsum。
假设 A 和 B 是两个长度相同的一维数组。例如A = np.arange(10)和B = np.arange(5, 15)。
的总和
A可以写成：
np.einsum('i->', A)
元素乘法
A * B可以写成：
np.einsum('i,i->i', A, B)
内积或点积，
np.inner(A, B)或np.dot(A, B)，可以写成：
np.einsum('i,i->', A, B) # or just use 'i,i'
外积
np.outer(A, B)可以写成：
np.einsum('i,j->ij', A, B)
对于二维数组，C和D，假设轴的长度兼容（两个长度相同或其中一个的长度为 1），这里有几个例子：
C（主对角线和）的迹np.trace(C)可以写成：
np.einsum('ii', C)
的元素乘法
C和 的转置D可以C * D.T写成：
np.einsum('ij,ji->ij', C, D)
将 的每个元素
C乘以数组D（形成四维数组），C[:, :, None, None] * D可以写成：
np.einsum('ij,kl->ijkl', C, D)
解决方案 2：
如果你直观地理解，那么掌握 的概念numpy.einsum()就很容易了。作为例子，让我们从涉及矩阵乘法的简单描述开始。
要使用numpy.einsum()，您只需将所谓的下标字符串作为参数传递，然后输入数组。
假设您有两个 2D 数组A和B，并且您想要进行矩阵乘法。因此，您可以这样做：
np.einsum("ij, jk -> ik", A, B)
这里下标字符串 ij对应于数组，A而下标字符串 jk对应于数组B。此外，这里要注意的最重要的一点是，每个下标字符串中的字符数必须与数组的维度匹配（即，二维数组为两个字符，三维数组为三个字符，依此类推）。如果你在下标字符串之间重复字符（在我们的例子中），那么这意味着你希望总和沿着这些维度进行。因此，它们的总和将减少（即，该维度将消失）。 j**ein
此符号后的下标字符串表示->结果数组的维度。如果将其留空，则所有内容将相加，并返回一个标量值作为结果。否则，结果数组将根据下标字符串具有维度。在我们的示例中，它将是ik。这是直观的，因为我们知道要使矩阵乘法起作用，数组中的列数A必须与数组中的行数匹配B，这就是这里发生的事情（即，我们通过重复下标字符串j中的字符来编码此知识）
np.einsum()这里还有一些示例，简洁地说明了在实现一些常见张量或nd 数组运算时的用途/功能。
输入
# a vector
In [197]: vec
Out[197]: array([0, 1, 2, 3])
# an array
In [198]: A
Out[198]:
array([[11, 12, 13, 14],
[21, 22, 23, 24],
[31, 32, 33, 34],
[41, 42, 43, 44]])
# another array
In [199]: B
Out[199]:
array([[1, 1, 1, 1],
[2, 2, 2, 2],
[3, 3, 3, 3],
[4, 4, 4, 4]])
1）矩阵乘法（类似np.matmul(arr1, arr2)）
In [200]: np.einsum("ij, jk -> ik", A, B)
Out[200]:
array([[130, 130, 130, 130],
[230, 230, 230, 230],
[330, 330, 330, 330],
[430, 430, 430, 430]])
2）沿主对角线提取元素（类似np.diag(arr)）
In [202]: np.einsum("ii -> i", A)
Out[202]: array([11, 22, 33, 44])
3）Hadamard 积（即两个数组的元素乘积）（类似arr1 * arr2）
In [203]: np.einsum("ij, ij -> ij", A, B)
Out[203]:
array([[ 11, 12, 13, 14],
[ 42, 44, 46, 48],
[ 93, 96, 99, 102],
[164, 168, 172, 176]])
4）逐元素平方（类似于np.square(arr)或arr ** 2）
In [210]: np.einsum("ij, ij -> ij", B, B)
Out[210]:
array([[ 1, 1, 1, 1],
[ 4, 4, 4, 4],
[ 9, 9, 9, 9],
[16, 16, 16, 16]])
5）迹（即主对角线元素之和）（类似于np.trace(arr)）
In [217]: np.einsum("ii -> ", A)
Out[217]: 110
6）矩阵转置（类似np.transpose(arr)）
In [221]: np.einsum("ij -> ji", A)
Out[221]:
array([[11, 21, 31, 41],
[12, 22, 32, 42],
[13, 23, 33, 43],
[14, 24, 34, 44]])
7) 向量外积（类似于np.outer(vec1, vec2)）
In [255]: np.einsum("i, j -> ij", vec, vec)
Out[255]:
array([[0, 0, 0, 0],
[0, 1, 2, 3],
[0, 2, 4, 6],
[0, 3, 6, 9]])
8) 向量内积（类似于np.inner(vec1, vec2)）
In [256]: np.einsum("i, i -> ", vec, vec)
Out[256]: 14
9）沿 0 轴求和（类似np.sum(arr, axis=0)）
In [260]: np.einsum("ij -> j", B)
Out[260]: array([10, 10, 10, 10])
10）沿轴 1 求和（类似np.sum(arr, axis=1)）
In [261]: np.einsum("ij -> i", B)
Out[261]: array([ 4, 8, 12, 16])
11）批量矩阵乘法
In [287]: BM = np.stack((A, B), axis=0)
In [288]: BM
Out[288]:
array([[[11, 12, 13, 14],
[21, 22, 23, 24],
[31, 32, 33, 34],
[41, 42, 43, 44]],
[[ 1, 1, 1, 1],
[ 2, 2, 2, 2],
[ 3, 3, 3, 3],
[ 4, 4, 4, 4]]])
In [289]: BM.shape
Out[289]: (2, 4, 4)
# batch matrix multiply using einsum
In [292]: BMM = np.einsum("bij, bjk -> bik", BM, BM)
In [293]: BMM
Out[293]:
array([[[1350, 1400, 1450, 1500],
[2390, 2480, 2570, 2660],
[3430, 3560, 3690, 3820],
[4470, 4640, 4810, 4980]],
[[ 10, 10, 10, 10],
[ 20, 20, 20, 20],
[ 30, 30, 30, 30],
[ 40, 40, 40, 40]]])
In [294]: BMM.shape
Out[294]: (2, 4, 4)
12）沿轴 2 求和（类似np.sum(arr, axis=2)）
In [330]: np.einsum("ijk -> ij", BM)
Out[330]:
array([[ 50, 90, 130, 170],
[ 4, 8, 12, 16]])
13）对数组中的所有元素求和（类似于np.sum(arr)）
In [335]: np.einsum("ijk -> ", BM)
Out[335]: 480
14）对多个轴求和（即边缘化）
（类似np.sum(arr, axis=(axis0, axis1, axis2, axis3, axis4, axis6, axis7))）
# 8D array
In [354]: R = np.random.standard_normal((3,5,4,6,8,2,7,9))
# marginalize out axis 5 (i.e. "n" here)
In [363]: esum = np.einsum("ijklmnop -> n", R)
# marginalize out axis 5 (i.e. sum over rest of the axes)
In [364]: nsum = np.sum(R, axis=(0,1,2,3,4,6,7))
In [365]: np.allclose(esum, nsum)
Out[365]: True
15）双点积（类似于np.sum(hadamard-product)参见3）
In [772]: A
Out[772]:
array([[1, 2, 3],
[4, 2, 2],
[2, 3, 4]])
In [773]: B
Out[773]:
array([[1, 4, 7],
[2, 5, 8],
[3, 6, 9]])
In [774]: np.einsum("ij, ij -> ", A, B)
Out[774]: 124
16）二维和三维数组乘法
当您求解线性方程组（ Ax = b）并想要验证结果时，这种乘法非常有用。
# inputs
In [115]: A = np.random.rand(3,3)
In [116]: b = np.random.rand(3, 4, 5)
# solve for x
In [117]: x = np.linalg.solve(A, b.reshape(b.shape[0], -1)).reshape(b.shape)
# 2D and 3D array multiplication :)
In [118]: Ax = np.einsum('ij, jkl', A, x)
# indeed the same!
In [119]: np.allclose(Ax, b)
Out[119]: True
相反，如果必须使用np.matmul()此验证，我们必须执行几个reshape操作才能达到相同的结果，例如：
# reshape 3D array `x` to 2D, perform matmul
# then reshape the resultant array to 3D
In [123]: Ax_matmul = np.matmul(A, x.reshape(x.shape[0], -1)).reshape(x.shape)
# indeed correct!
In [124]: np.allclose(Ax, Ax_matmul)
Out[124]: True
奖励：在这里阅读更多数学知识：Einstein-Summation，当然还有这里：Tensor-Notation
解决方案 3：
当阅读 einsum 方程时，我发现最有帮助的是能够在脑海中将它们归结为它们的命令式版本。
让我们从以下（强加的）声明开始：
C = np.einsum('bhwi,bhwj->bij', A, B)
首先查看标点符号，我们发现箭头前有两个 4 个字母的逗号分隔的斑点 -bhwi和，箭头后有一个 3 个字母的斑点。因此，该方程从两个 4 阶张量输入生成一个 3 阶张量结果。bhwj`bij`
现在，让每个 blob 中的每个字母成为范围变量的名称。字母在 blob 中出现的位置是它在该张量中范围的轴的索引。因此，产生 C 的每个元素的命令式求和必须以三个嵌套的 for 循环开始，每个 C 的索引一个。
for b in range(...):
for i in range(...):
for j in range(...):
# the variables b, i and j index C in the order of their appearance in the equation
C[b, i, j] = ...
因此，本质上，您for对 C 的每个输出索引都有一个循环。我们暂时不确定范围。
接下来我们看看左侧 - 那里是否有任何未出现在右侧的范围变量？ 在我们的例子中 - 是的，h并且w。 为每个这样的变量添加一个内部嵌套for循环：
for b in range(...):
for i in range(...):
for j in range(...):
C[b, i, j] = 0
for h in range(...):
for w in range(...):
...
在最内层循环中，我们现在已经定义了所有索引，因此我们可以写出实际的总和，并且翻译就完成了：
# three nested for-loops that index the elements of C
for b in range(...):
for i in range(...):
for j in range(...):
# prepare to sum
C[b, i, j] = 0
# two nested for-loops for the two indexes that don't appear on the right-hand side
for h in range(...):
for w in range(...):
# Sum! Compare the statement below with the original einsum formula
# 'bhwi,bhwj->bij'
C[b, i, j] += A[b, h, w, i] * B[b, h, w, j]
如果您到目前为止能够理解代码，那么恭喜您！这就是您能够理解 einsum 方程式所需的全部内容。特别注意原始 einsum 公式如何映射到上面代码片段中的最终求和语句。for 循环和范围边界只是虚词，而最终语句才是您真正需要了解发生了什么的全部内容。
为了完整起见，让我们看看如何确定每个范围变量的范围。好吧，每个变量的范围只是它索引的维度的长度。显然，如果一个变量在一个或多个张量中索引多个维度，那么每个维度的长度必须相等。以下是上面带有完整范围的代码：
# C's shape is determined by the shapes of the inputs
# b indexes both A and B, so its range can come from either A.shape or B.shape
# i indexes only A, so its range can only come from A.shape, the same is true for j and B
assert A.shape[0] == B.shape[0]
assert A.shape[1] == B.shape[1]
assert A.shape[2] == B.shape[2]
C = np.zeros((A.shape[0], A.shape[3], B.shape[3]))
for b in range(A.shape[0]): # b indexes both A and B, or B.shape[0], which must be the same
for i in range(A.shape[3]):
for j in range(B.shape[3]):
# h and w can come from either A or B
for h in range(A.shape[1]):
for w in range(A.shape[2]):
C[b, i, j] += A[b, h, w, i] * B[b, h, w, j]
解决方案 4：
另一种观点np.einsum
这里的大多数答案都是通过例子来解释的，我想我会给出一个额外的观点。
您可以将其视为einsum广义矩阵求和运算符。给定的字符串包含下标，它们是代表轴的标签。我喜欢将其视为您的操作定义。下标提供了两个明显的约束：
每个输入数组的轴数，
输入之间的轴大小相等。
让我们以最初的例子为例：np.einsum('ij,jk->ki', A, B)。这里的约束1.转化为 A.ndim == 2和B.ndim == 2，以及2.转化为A.shape[1] == B.shape[0]。
正如您稍后会看到的，还有其他限制。例如：
输出下标中的标签不得出现超过一次。
输出下标中的标签必须出现在输入下标中。
当查看时ij,jk->ki，您可以将其视为：
输入数组中的哪些组件将会对
[k, i]输出数组的组件产生贡献。
下标包含输出数组每个组件的操作的精确定义。
我们将坚持操作，以及和ij,jk->ki的以下定义：A`B`
>>> A = np.array([[1,4,1,7], [8,1,2,2], [7,4,3,4]])
>>> A.shape
(3, 4)
>>> B = np.array([[2,5], [0,1], [5,7], [9,2]])
>>> B.shape
(4, 2)
输出Z的形状为 ，(B.shape[1], A.shape[0])可以简单地按以下方式构造。从 的空白数组开始Z：
Z = np.zeros((B.shape[1], A.shape[0]))
for i in range(A.shape[0]):
for j in range(A.shape[1]):
for k range(B.shape[0]):
Z[k, i] += A[i, j]*B[j, k] # ki <- ij*jk
np.einsum是关于在输出数组中累积贡献。每对(A[i,j], B[j,k])都对每个组件有贡献Z[k, i]。
您可能已经注意到，它看起来与计算一般矩阵乘法的方式极为相似......
最低限度的实施
这是 Python 中的最小实现np.einsum。这应该有助于理解幕后真正发生的事情。
随着我们的讨论，我将继续参考前面的例子。定义inputs为[A, B]。
np.einsum实际上可以接受两个以上的输入。下面，我们将重点介绍一般情况：n 个输入和n 个输入下标。主要目标是找到迭代域，即所有范围的笛卡尔积。
我们不能依赖手动编写for循环，因为我们不知道会有多少个循环。主要思想是这样的：我们需要找到所有唯一标签（我将使用key和keys来引用它们），找到相应的数组形状，然后为每个标签创建范围，并计算范围的乘积itertools.product以获得研究域。
| 指数 | keys | 约束 | sizes | ranges |
|---|---|---|---|---|
| 1 | 'i' | A.shape[0] | 3 | range(0, 3) |
| 2 | 'j' | A.shape[1] == B.shape[0] | 4 | range(0, 4) |
| 0 | 'k' | B.shape[1] | 2 | range(0, 2) |
研究的范围是笛卡尔积：range(0, 2) x range(0, 3) x range(0, 4)。
下标处理：
>>> expr = 'ij,jk->ki'
>>> qry_expr, res_expr = expr.split('->')
>>> inputs_expr = qry_expr.split(',')
>>> inputs_expr, res_expr
(['ij', 'jk'], 'ki')
查找输入下标中的唯一键（标签）：
>>> keys = set([(key, size) for keys, input in zip(inputs_expr, inputs)
for key, size in list(zip(keys, input.shape))])
{('i', 3), ('j', 4), ('k', 2)}
我们应该检查约束（以及输出下标）！使用set不是一个好主意，但它对于本例的目的来说是可行的。
获取相关大小（用于初始化输出数组）并构造范围（用于创建我们的迭代域）：
>>> sizes = dict(keys)
{'i': 3, 'j': 4, 'k': 2}
>>> ranges = [range(size) for _, size in keys]
[range(0, 2), range(0, 3), range(0, 4)]
我们需要一个包含键（标签）的列表：
>>> to_key = sizes.keys()
['k', 'i', 'j']
range计算s的笛卡尔积
>>> domain = product(*ranges)
注意：[itertools.product][1]返回一个随着时间推移而被消耗的迭代器。
将输出张量初始化为：
>>> res = np.zeros([sizes[key] for key in res_expr])
我们将循环
domain：
>>> for indices in domain:
... pass
对于每次迭代，indices将包含每个轴上的值。在我们的示例中，这将提供i、j和k作为元组：(k, i, j)。对于每个输入（A和B），我们需要确定要获取哪个组件。没错，就是A[i, j]和B[j, k]！但是，从字面上讲，我们没有变量i、j和k。
我们可以indices使用zip 来创建每个键（标签）和其当前值to_key之间的映射：
>>> vals = dict(zip(to_key, indices))
为了获取输出数组的坐标，我们使用vals并循环遍历键：[vals[key] for key in res_expr]。但是，要使用这些键来索引输出数组，我们需要用tuple和zip来包装它以分隔每个轴上的索引：
>>> res_ind = tuple(zip([vals[key] for key in res_expr]))
输入索引相同（尽管可以有多个）：
>>> inputs_ind = [tuple(zip([vals[key] for key in expr])) for expr in inputs_expr]
我们将使用 a
itertools.reduce来计算所有贡献成分的乘积：
>>> def reduce_mult(L):
... return reduce(lambda x, y: x*y, L)
总体而言，域上的循环如下所示：
>>> for indices in domain:
... vals = {k: v for v, k in zip(indices, to_key)}
... res_ind = tuple(zip([vals[key] for key in res_expr]))
... inputs_ind = [tuple(zip([vals[key] for key in expr]))
... for expr in inputs_expr]
...
... res[res_ind] += reduce_mult([M[i] for M, i in zip(inputs, inputs_ind)])
>>> res
array([[70., 44., 65.],
[30., 59., 68.]])
这与np.einsum('ij,jk->ki', A, B)回报非常接近！
解决方案 5：
我发现《NumPy：行业秘诀（第二部分）》很有启发
我们使用 -> 来表示输出数组的顺序。因此，将“ij, i->j”视为具有左侧 (LHS) 和右侧 (RHS)。LHS 上的任何重复标签都会按元素计算乘积，然后求和。通过更改 RHS（输出）侧的标签，我们可以定义我们想要相对于输入数组进行的轴，即沿轴 0、1 等求和。
import numpy as np
>>> a
array([[1, 1, 1],
[2, 2, 2],
[3, 3, 3]])
>>> b
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> d = np.einsum('ij, jk->ki', a, b)
注意有三个轴，i、j、k，并且 j 重复（在左侧）。 i,j表示 的行和列a。表示j,k的列b。
为了计算乘积并对齐j轴，我们需要添加一个轴a。（b将沿着第一个轴广播？）
a[i, j, k]
b[j, k]
>>> c = a[:,:,np.newaxis] * b
>>> c
array([[[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8]],
[[ 0, 2, 4],
[ 6, 8, 10],
[12, 14, 16]],
[[ 0, 3, 6],
[ 9, 12, 15],
[18, 21, 24]]])
j右侧不存在，因此我们将其相加，j得到 3x3x3 数组的第二个轴
>>> c = c.sum(1)
>>> c
array([[ 9, 12, 15],
[18, 24, 30],
[27, 36, 45]])
最后，索引（按字母顺序）在右侧反转，因此我们进行转置。
>>> c.T
array([[ 9, 18, 27],
[12, 24, 36],
[15, 30, 45]])
>>> np.einsum('ij, jk->ki', a, b)
array([[ 9, 18, 27],
[12, 24, 36],
[15, 30, 45]])
>>>
解决方案 6：
让我们制作两个具有不同但兼容维度的数组，以突出它们的相互作用
In [43]: A=np.arange(6).reshape(2,3)
Out[43]:
array([[0, 1, 2],
[3, 4, 5]])
In [44]: B=np.arange(12).reshape(3,4)
Out[44]:
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
您的计算是，将 (2,3) 与 (3,4) 的“点”（乘积之和）相加，得到 (4,2) 数组。 i是 的第一个维数A， 的最后一个维数C；k的最后一个维数B， 的第一个维数C。 j被求和“消耗”。
In [45]: C=np.einsum('ij,jk->ki',A,B)
Out[45]:
array([[20, 56],
[23, 68],
[26, 80],
[29, 92]])
np.dot(A,B).T这与转置后的最终输出相同。
要查看 发生的更多信息j，请将C下标更改为ijk：
In [46]: np.einsum('ij,jk->ijk',A,B)
Out[46]:
array([[[ 0, 0, 0, 0],
[ 4, 5, 6, 7],
[16, 18, 20, 22]],
[[ 0, 3, 6, 9],
[16, 20, 24, 28],
[40, 45, 50, 55]]])
这也可以通过以下方式生成：
A[:,:,None]*B[None,:,:]
也就是说，k在 的末尾添加一个维度A，i在 的前面添加一个维度B，得到一个 (2,3,4) 数组。
0 + 4 + 16 = 20，，9 + 28 + 55 = 92等；求和j并转置得到先前的结果：
np.sum(A[:,:,None] * B[None,:,:], axis=1).T
# C[k,i] = sum(j) A[i,j (,k) ] * B[(i,) j,k]
解决方案 7：
一旦熟悉了虚拟索引（公共或重复索引）以及爱因斯坦求和（einsum）中沿虚拟索引的求和，输出->整形就很容易了。因此，重点关注：
虚拟索引，
j是np.einsum("ij,jk->ki", a, b)沿虚拟索引求和
j
虚拟索引
对于，无论是否有共同的指标，einsum("...", a, b)元素乘法总是发生在矩阵a和之间。我们可以得到在下标 中没有共同指标的。b`einsum('xy,wz', a, b)`'xy,wz'
如果有一个公共指标，如j，那么在爱因斯坦求和中"ij,jk->ki"它被称为虚拟指标。
爱因斯坦求和
被求和的指标是求和指标，在本例中为“i”。它也被称为虚拟指标，因为任何符号都可以替换“i”，而不会改变表达式的含义，只要它不与同一术语中的指标符号冲突。
沿虚拟索引求和
对于图中绿色矩形np.einsum("ij,j", a, b)的，是虚拟索引。元素乘法沿轴相加为。j`ai b[j]jΣ ( ai b[j] )`
它是每个的点积 。这里具体使用和避免，因为它不是严格意义上的数学点积运算。np.inner(a[i], b)`inp.inner()np.dot`
爱因斯坦求和约定：简介
只要满足规则（详情请参阅youtube），虚拟索引可以出现在任何位置。
i对于中的虚拟索引np.einsum(“ik,il", a, b)，它是矩阵a和的行索引b，因此从a 和 中b提取一列来生成点积s。
输出形式
因为求和是沿着虚拟索引进行的，所以虚拟索引在结果矩阵中消失，因此i从“ik,il"被删除并形成形状(k,l)。我们可以通过带有标识符的输出下标标签np.einsum("... -> <shape>")来指定输出形式。->
有关详细信息，请参阅numpy.einsum中的显式模式。
在显式模式下，可以通过指定输出下标标签来直接控制输出。这需要标识符
‘->’以及输出下标标签列表。此功能增加了函数的灵活性，因为可以在需要时禁用或强制求和。调用np.einsum('i->', a)类似于np.sum(a, axis=-1)，并且np.einsum('ii->i', a)类似于np.diag(a)。不同之处在于 einsum 默认不允许广播。此外，np.einsum('ij,jh->ih', a, b)直接指定输出下标标签的顺序，因此返回矩阵乘法，这与上面的隐式模式示例不同。
无虚拟索引
einsum 中没有虚拟索引的一个例子。
一个术语（下标索引，例如
"ij"）选择每个数组中的一个元素。每个左侧元素都应用于右侧元素进行逐元素乘法（因此总是发生乘法）。
a具有形状 (2,3)，其中每个元素应用于b形状 (2,2)。因此，它创建了一个形状(2,3,2,2)为无求和的矩阵(i,j)，因为(k.l)都是自由指标。
# --------------------------------------------------------------------------------
# For np.einsum("ij,kl", a, b)
# 1-1: Term "ij" or (i,j), two free indices, selects selects an element a[i][j].
# 1-2: Term "kl" or (k,l), two free indices, selects selects an element b[k][l].
# 2: Each a[i][j] is applied on b[k][l] for element-wise multiplication a[i][j] * b[k,l]
# --------------------------------------------------------------------------------
# for (i,j) in a:
# for(k,l) in b:
# a[i][j] * b[k][l]
np.einsum("ij,kl", a, b)
array([[[[ 0, 0],
[ 0, 0]],
[[10, 11],
[12, 13]],
[[20, 22],
[24, 26]]],
[[[30, 33],
[36, 39]],
[[40, 44],
[48, 52]],
[[50, 55],
[60, 65]]]])
示例
矩阵 A 行和矩阵 B 列的点积
A = np.matrix('0 1 2; 3 4 5')
B = np.matrix('0 -3; -1 -4; -2 -5');
np.einsum('ij,ji->i', A, B)
# Same with
np.diagonal(np.matmul(A,B))
(A*B).diagonal()
---
[ -5 -50]
[ -5 -50]
[[ -5 -50]]
解决方案 8：
我认为最简单的例子在tensorflow 文档中
将方程式转换为 einsum 符号有四个步骤。让我们以此方程式为例C[i,k] = sum_j A[i,j] * B[j,k]
首先我们删除变量名。我们得到
ik = sum_j ij * jk我们删除该项
sum_j，因为它是隐式的。我们得到ik = ij * jk我们
*用替换,。我们得到ik = ij, jk输出在 RHS 上，并用
->符号分隔。我们得到ij, jk -> ik
einsum 解释器只是反向运行这 4 个步骤。结果中缺失的所有索引都会被加起来。
以下是文档中的更多示例
# Matrix multiplication
einsum('ij,jk->ik', m0, m1) # output[i,k] = sum_j m0[i,j] * m1[j, k]
# Dot product
einsum('i,i->', u, v) # output = sum_i u[i]*v[i]
# Outer product
einsum('i,j->ij', u, v) # output[i,j] = u[i]*v[j]
# Transpose
einsum('ij->ji', m) # output[j,i] = m[i,j]
# Trace
einsum('ii', m) # output[j,i] = trace(m) = sum_i m[i, i]
# Batch matrix multiplication
einsum('aij,ajk->aik', s, t) # out[a,i,k] = sum_j s[a,i,j] * t[a, j, k]
