【NumPy】爱因斯坦求和约定

标记法

NumPy方法名称einsum，其全称是Einstein summation convention（爱因斯坦求和约定），又称为爱因斯坦标记法（Einstein notation）。在处理关于坐标的方程式时，这个方法非常有用。这是由大名鼎鼎的物理学家阿尔伯特·爱因斯坦（Albert Einstein）于1916年提出的。
这种约定，简单来说，就是省去了求和式中的求和符号。我们来看下面的点乘公式。
$$s=\sum _i{v}_i{w}_i$$
这个公式表述的是，两个向量对应的元素相乘后求和。这种写法于我们而言，已然熟稔于心。可爱因斯坦偏偏觉得这种数学符号太过于烦琐，他觉得那个求和符号纯属多余，于是他发明了另外一种写法。
$$s=v_i{w}^i$$
请特别注意，在爱因斯坦的标记体系里，下标表示行向量中的元素：
KaTeX parse error: Undefined control sequence: \cdts at position 10: [v_1,v_2,\̲c̲d̲t̲s̲,v_k]
但上标并不表示指数，而表示列向量中的元素：
$$\left[\begin{array}{c}{w}^1\\ w^2\\ \cdots \\ w^k\end{array}\right]$$

1. 内积：
   $$u\cdot v=u_j{v}^i$$
2. 向量乘以矩阵：矩阵$A$和向量$v$的乘积向量$u$
   $$u^i=A_j^i{v}^j$$
3. 矩阵乘法：矩阵$A_{i\ast j}$和矩阵$B_{j\ast k}$，二者乘法
   $$C_k^i=A_j^i{B}_k^j$$
4. 矩阵的迹：对于一个矩阵A，如果矩阵的上标与下标相同，则可以得到这个矩阵的迹t。
   $$t=A_i^i$$
5. 外积：M维向量a和N维余向量b的外积是一个M×N的矩阵A。
   $$A=ab$$

如果采用爱因斯坦标记法，上述公式可以表示如下。
$$A_j^i=a^i{b}_j$$
由于i和j代表两个不同的标号，外积不会除去这两个标号，于是这两个标号变成了新矩阵A的标号。

NumPy中的einsum()方法

格式：

result = np.einsum("{}, {}, {}->{}", arg0, arg1, arg2)

einsum()方法，在理论上，可以支持任意多的参数。该方法的第一个参数—格式字符串，至关重要，它直接决定整个方法实现的功能。格式字符串的规定如下。

• 不同的输入变量格式字符之间要用逗号分隔开，且输入格式字符的数量要与参与运算的变量数量相匹配。
• 输入格式字符和输出格式字符要用箭头分隔开。
• 输入格式字符和输出格式字符都是一系列常见字符（如ASCII码）。
• 格式字符串中的字符数（字符串长度）和张量的维数相对应。

np.einsum()的用法

import numpy as np

arr = np.arange(10)
print("\narr: \n{}".format(arr))
sum1 = np.einsum('i->', arr)
print("\nsum1: {}".format(sum1))
======================================
arr: 
[0 1 2 3 4 5 6 7 8 9]

sum1: 45

格式字符串'i->'表示的是一个一维张量（箭头前面的字符长度为1），它的维度消失了（箭头后面的字符为空，其实是变成了一个标量），怎么才能达到这个效果呢？再结合einsum()方法的关键词“sum”可知，原来是通过“求和”达成的“约减”。

按行求和 按列求和

如果按照上述逻辑，巧妙利用einsum()中的格式字符，我们可以很容易实现“按行求和”或“按列求和”的操作。

arr2 = np.arange(20).reshape(4, 5)
print("\narr2:\n{}".format(arr2))
sum_col = np.einsum('ij->j', arr2)
print("\nsum_col:\n{}".format(sum_col))
============================================
arr2:
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]]

sum_col:
[30 34 38 42 46]

上述代码的核心依然是einsum()方法的格式字符串，在’ij->j’中，箭头之前的输入格式字符串由两个字符构成，说明这是一个二维张量，这和arr的特征是一致的。它已经变成一个4行5列的矩阵。按照字符出现的顺序，第一个字符“i”表示行，第二个字符“j”表示列，如前面的介绍，箭头之后的格式字符“j”表示输出变量的形态，可以发现两个变化：第一，字符数量减少为1个，这说明输出变量被“降维”了；第二，原来的字符“j”被保留，它原来的位置表示列，这说明“行”这个维度被消灭了。于是，放在一起，它的效果就是通过求和的方式达到维度约减。用学术点的话来说，上述操作实现了按列求和.

sum_row = np.einsum("ab->a", arr2)
print("\nsum_row:\n{}".format(sum_row))
============================================
sum_row:
[10 35 60 85]

在上面的代码中，我们故意把格式字符串变化为"ab->a"，就是想告诉你，在格式字符串中用什么字符并不重要，重要的是它们的长度，以及它们是否在箭头之后出现。上述格式字符串完全等价于"ij->i"。类似前面的分析，In [10]处实现的功能依然是通过求和进而降维。所不同的是，第一个字符a（代表行）保留，而第二个字符b（代表列）消失，这个格式表明，这是一个按行求和操作.

前面仅仅提到了利用einsum()求和，实际上，该方法的功能远不止于此。我们还可以利用它实现矩阵乘法，因为矩阵乘法涉及求和。

A = np.array([[1, 1, 1], [2, 2, 2], [5, 5, 5]])
B = np.array([[0, 1, 0], [1, 1, 0], [1, 1, 1]])
result = A @ B
print("\nresult:\n{}".format(result))
result2 = np.einsum('ij, jk->ik', A, B)
print("\nresult2:\n{}".format(result2))
===============================================
result:
[[ 2  3  1]
 [ 4  6  2]
 [10 15  5]]

result2:
[[ 2  3  1]
 [ 4  6  2]
 [10 15  5]]

利用einsum()做矩阵乘法，关键还是在于格式字符串。在上述代码的In [17]处，数组A对应的格式字符串为“ij”，字符长度为2，说明它是二维数组。类似地，二维数组B对应的格式字符串为“jk”，它们有相同的字符“j”。根据“j”出现的顺序，它表明第一个数组A的列和第二个数组B的行是相同的。但通过计算的结果发现，格式字符串中的箭头之后少了“j”的身影，这表明通过求和计算将“j”降维了，只有矩阵乘法才有类似的功效.

向量乘积

我们想求得两个向量（向量a和向量a）的中对应元素的乘积（Element-wisemultiplication），即对向量实施点乘操作，根据前面所讲，可以有两种写法：a*b或np.multiply(a,b)。学习了einsum()之后，还可以有第三种方法

a = np.array([[1, 2], [3, 4]])
b = np.ones(shape=(2, 2))
print("\narray:\n{}".format(np.einsum('ij,ij->ij', a, b)))
====================================================
array:
[[1. 2.]
 [3. 4.]]

向量内积

如果向量a和向量b是一维向量，则可写成np.einsum(‘i,i->i’,a,b)。
求向量a和向量b的内积时，可用np.inner(a, b)实现，也可以用einsum()轻易实现。

print(np.einsum('ij,ij->', a, b))
======================
10.0

矩阵的对角线

AA = np.array([[11, 12, 13, 14],
               [21, 22, 23, 24],
               [31, 32, 33, 34],
               [41, 42, 43, 44]])
print("\narray:\n{}".format(np.einsum('ii->i', AA)))
===========================================
array:
[11 22 33 44]

矩阵的迹

print(np.einsum('ii->', AA))
===================================
110

矩阵的转置

print(np.einsum('ij->ji', AA))
========================================
[[11 21 31 41]
 [12 22 32 42]
 [13 23 33 43]
 [14 24 34 44]]

完整代码

import numpy as np

arr = np.arange(10)
print("\narr: \n{}".format(arr))
sum1 = np.einsum('i->', arr)
print("\nsum1: {}".format(sum1))

arr2 = np.arange(20).reshape(4, 5)
print("\narr2:\n{}".format(arr2))
sum_col = np.einsum('ij->j', arr2)
print("\nsum_col:\n{}".format(sum_col))

sum_row = np.einsum("ab->a", arr2)
print("\nsum_row:\n{}".format(sum_row))

A = np.array([[1, 1, 1], [2, 2, 2], [5, 5, 5]])
B = np.array([[0, 1, 0], [1, 1, 0], [1, 1, 1]])
result = A @ B
print("\nresult:\n{}".format(result))
result2 = np.einsum('ij, jk->ik', A, B)
print("\nresult2:\n{}".format(result2))

a = np.array([[1, 2], [3, 4]])
b = np.ones(shape=(2, 2))
print("\narray:\n{}".format(np.einsum('ij,ij->ij', a, b)))

print(np.einsum('ij,ij->', a, b))

AA = np.array([[11, 12, 13, 14],
               [21, 22, 23, 24],
               [31, 32, 33, 34],
               [41, 42, 43, 44]])
print("\narray:\n{}".format(np.einsum('ii->i', AA)))

print(np.einsum('ii->', AA))
print(np.einsum('ij->ji', AA))