我发现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，所以我们对3x3x3数组的第二个轴j求和

>>> 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]])
>>>

让我们创建两个数组，它们具有不同但兼容的尺寸，以突出它们的相互作用

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的第一个dim， C的最后一个dim;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,:,:]

也就是说，在A的末尾添加一个k维度，在B的前面添加一个i维度，从而得到一个(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]

(注意:这个答案是基于我之前写的关于einsum的简短博客。)

einsum做什么?

假设我们有两个多维数组A和B。现在让我们假设我们想……

• 乘 A与B以特定的方式创建新的产品数组;然后也许
• 总和这个新数组沿着特定的轴;然后也许
• 转置新数组的轴以特定的顺序。

相比multiply、sum和transpose等NumPy函数的组合，einsum很有可能帮助我们更快、更有效地执行此操作。

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(我们希望最终得到一个1D数组)。通过省略标签，我们将沿着这个轴告诉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被保留用于输出，因此我们得到一个2D数组。

将这个结果与标签j为不求和的数组进行比较可能会更清楚。下面，在左边你可以看到写入np.einsum('ij,jk->ijk', A, B)(即我们保留标签j)后的3D数组:

求和轴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)
  

对于2D数组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(生成一个4D数组)，C[:, :, None, None] * D，可以写成:

  np.einsum('ij,kl->ijkl', C, D)
  

如果你能直观地理解< >强numpy.einsum() < / >强，就很容易理解它。作为一个例子，让我们从一个涉及矩阵乘法的简单描述开始。

要使用< >强numpy.einsum() < / >强，你所要做的就是传递所谓的下标字符串作为参数，后面跟着你的输入数组。

假设你有两个2D数组，A和B，你想做矩阵乘法。所以，你会:

np.einsum("ij, jk -> ik", A, B)

这里的下标字符串 ij对应数组A，而下标字符串 A0对应数组A1。此外，这里要注意的最重要的事情是，每个下标字符串 A4中的A2与数组的尺寸匹配(即，2D数组为两个字符，3D数组为三个字符，等等)。如果你在A5(在我们的例子中是A6)之间重复字符，那么这意味着你希望einA7在这些维度上发生。因此，它们将被求和-约简(即，该维度将是A8)。

下面是一些例子，简洁地说明了< >强np.einsum() < / >强在实现一些常见的张量或nd-array操作时的使用/功能。

输入

# 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) Trace(即主对角元素之和)(类似于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

# 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) cf. 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) 2D和3D数组乘法

这样的乘法在求解线性方程组(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

奖金:在这里阅读更多数学:爱因斯坦求和和绝对在这里:Tensor-Notation

让我们从下面的陈述开始:

C = np.einsum('bhwi,bhwj->bij', A, B)

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] = ...

因此，本质上，对于c的每个输出索引，都有一个for循环。我们暂时不确定范围。

接下来我们看左边——有没有任何范围变量不出现在右手边?在我们的例子中，是的，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]

我认为最简单的例子是tensorflow文档

将方程转换为einsum符号有四个步骤。让我们以这个方程为例C[i,k] = sum_j A[i,j] * B[j,k]

1. 首先，我们删除变量名。我们得到ik = sum_j ij * jk
2. 我们删除sum_j项，因为它是隐式的。我们得到ik = ij * jk
3. 我们将*替换为,。我们得到ik = ij, jk
4. 输出在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]

np.einsum的另一个视图

这里的大多数答案都是通过例子来解释的，我想我会给出一个额外的观点。

你可以看到einsum是一个广义矩阵求和运算符。给定的字符串包含下标，下标是表示轴的标签。我喜欢把它看作你的操作定义。下标提供了两个明显的约束条件:

1. 每个输入数组的轴数，

2. 输入之间的轴大小相等。

让我们以最初的例子为例:np.einsum('ij,jk->ki', A, B)。在这里，约束<强> 1。< / >强转换为A.ndim == 2和B.ndim == 2，而< >强2。< / >强转换为A.shape[1] == B.shape[0]。

稍后您将看到，还有其他限制条件。例如:

3. 标签在输出下标中不能出现超过一次。

4. 输出下标中的标签必须出现在输入下标中。

当观察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创建范围的乘积计算来获得研究的域。

研究的领域是笛卡尔积:range(0, 2) x range(0, 3) x range(0, 4)。

1. < p >下标处理:

   >>> expr = 'ij,jk->ki'
   >>> qry_expr, res_expr = expr.split('->')
   >>> inputs_expr = qry_expr.split(',')
   >>> inputs_expr, res_expr
   (['ij', 'jk'], 'ki')
   

2. 找到输入下标中的唯一键(标签):

   >>> 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是一个坏主意，但它将适用于本例的目的。

3. 获取相关的大小(用于初始化输出数组)并构造范围(用于创建迭代域):

   >>> 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)]
   

4. 我们需要一个包含键的列表 (标签):

   >>> to_key = sizes.keys()
   ['k', 'i', 'j']
   

5. 计算__abc0的笛卡尔积

   >>> domain = product(*ranges)
   

   请注意: [itertools.product][1]返回一个迭代器，该迭代器随着时间的推移获得消耗。

6. 初始化输出张量为:

   >>> res = np.zeros([sizes[key] for key in res_expr])
   

7. 我们将循环domain:

   >>> for indices in domain:
   ...     pass
   

   对于每次迭代，indices将包含每个轴上的值。在我们的例子中，这将提供i， j和k作为i2: (k, i, j)。对于每个输入(A和B)，我们需要确定要获取的组件。那是A[i, j]和B[j, k]，对!然而，从字面上讲，我们没有变量i， j和k。

   我们可以用to_key压缩indices来创建每个键(标签)与其当前值之间的映射:

   >>> 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]
   

8. 我们将使用itertools.reduce来计算所有贡献组件的乘积:

   >>> def reduce_mult(L):
   ...     return reduce(lambda x, y: x*y, L)
   

9. 总的来说，在域上的循环看起来像:

   >>> 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)返回的结果非常接近!

一旦熟悉了虚拟索引(公共索引或重复索引)和爱因斯坦总结 (einsum)中虚拟索引的和，输出->的整形就很容易了。因此，重点关注:

1. 虚拟索引，即np.einsum("ij,jk->ki", a, b)中的公共索引j
2. 沿着虚索引j求和

哑指标

对于einsum("...", a, b)，元素乘总是发生在矩阵a和b之间，无论是否有公共索引。我们可以有einsum('xy,wz', a, b)，它在下标'xy,wz'中没有公共索引。

如果存在一个公共索引，如"ij,jk->ki"中的j，那么它在爱因斯坦求和中被称为哑指标。

• 爱因斯坦求和

被求和的索引是求和索引，在本例中是&;i&;。它也被称为虚拟索引，因为任何符号都可以取代“;i"不改变表达式的含义，前提是它不与同一项中的索引符号冲突。

沿虚指标求和

对于图中绿色矩形的np.einsum("ij,j", a, b)， j是虚拟索引。元素的乘法a[i][j] * b[j]沿着j轴求和为Σ ( a[i][j] * b[j] )。

每个i对应一个点积 np.inner(a[i], b)。这里具体到np.inner()，避免使用np.dot，因为它不是严格的数学点积操作。

• 爱因斯坦求和约定:一个介绍

虚拟索引可以出现在任何地方，只要规则满足(详情请参阅youtube)。

输出形式

因为求和是沿着哑指标进行的，所以虚拟索引在结果矩阵中消失，因此“ik,il"中的i被删除并形成形状(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中没有虚索引的例子。

1. 一个term(下标indexes，例如"ij")在每个数组中选择一个元素。
2. 左边的每个元素都应用到右边的元素上，用于按元素进行乘法(因此总是会发生乘法)。

a具有shape(2,3)，其中每个元素应用于shape(2,2)的b。因此，它创建了一个形状为(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]]