np。einsum工作吗?

给定数组A和B,它们的矩阵乘法和转置是使用(A @ B).T计算的,或者等效地,使用:

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

当前回答

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

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

也就是说,在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]

其他回答

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

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

也就是说,在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]

我认为最简单的例子是在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]

我发现NumPy: The tricks of The trade (Part II)很有教育意义

我们使用->表示输出数组的顺序。把'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求和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]])
>>>

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

虚索引,np中的公共索引j。Einsum ("ij,jk->ki", a, b) 沿虚指标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])。

它是np的点积。对于每个i, inner(a[i], b)是特定的np.inner()而避免np。因为它不是严格意义上的点积运算。

爱因斯坦求和惯例:简介

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

对于np中的虚索引i。Einsum (" ik,il", a, b),它是矩阵a和b的行索引,因此提取a的列和b的列来生成点积。

输出形式

由于求和是沿着假指标进行的,所以假指标在结果矩阵中消失,因此i从“ik,il”中去掉,形成形状(k,l)。我们可以告诉np.einsum(“…-> <shape>")通过带有->标识符的输出下标标签指定输出表单。

参见numpy中的显式模式。详情请参阅Einsum。

In explicit mode the output can be directly controlled by specifying output subscript labels. This requires the identifier ‘->’ as well as the list of output subscript labels. This feature increases the flexibility of the function since summing can be disabled or forced when required. The call np.einsum('i->', a) is like np.sum(a, axis=-1), and np.einsum('ii->i', a) is like np.diag(a). The difference is that einsum does not allow broadcasting by default. Additionally np.einsum('ij,jh->ih', a, b) directly specifies the order of the output subscript labels and therefore returns matrix multiplication, unlike the example above in implicit mode.

没有虚索引

一个在einsum中没有虚索引的例子。

一个术语(下标索引,例如:"ij")选择每个数组中的一个元素。 左边的每个元素都应用到右边的元素上,用于按元素进行乘法(因此总是会发生乘法)。

A具有形状(2,3),其中每个元素都应用于形状(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]]

在阅读einsum方程时,我发现能够这样做最有帮助 在心里把它们归结为它们的命令版本。

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

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

首先看标点符号,我们看到箭头前面有两个4个逗号分隔的斑点——bhwi和bhwj, 后面还有一个3个字母的bij。因此,该方程由两个秩4张量输入产生一个秩3张量结果。

现在,让每个blob中的每个字母都是一个范围变量的名称。字母出现在圆点中的位置 是它在这个张量中所覆盖的轴的下标。 因此,生成C的每个元素的命令式求和必须从三个嵌套的for循环开始,每个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] = ...

所以,本质上,对于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]