Python的Numpy提供了很多高效的科学计算函数,einsum便是其中一个。einsum函数非常灵活,也非常高效,运行时占用的存储空间很小。由于其强大的表现力和智能循环,它在速度和内存效率方面通常可以超越我们常见的array函数。但缺点是,可能需要一段时间才能理解符号,有时需要尝试才能将其正确的应用于棘手的问题。

1. 爱因斯坦求和约定

在线性代数中,我们最多涉及的是二阶及以下的张量.在这种情况下,纸面上可以很方便地写出低阶张量的矩阵形式,高阶的张量,它们的坐标就没法用矩阵表示.我们当然可以把矩阵拓展为立体阵等概念,但随着阶数上升,这种表示法的复杂程度几何级增加;我们也可以使用张量词条中所提过的向量矩阵的方法,比起立体阵要清楚一些,但套娃式的表达方式也对理解一个张量的性质造成了障碍.

爱因斯坦求和约定正是为了简洁地表达高阶张量的坐标运算而存在的.

假设 矩阵大小分别是 2*3 和 3*2,矩阵乘法的定义如下:

爱因斯坦求和是一种对求和公式简洁高效的记法,其原则是当变量下标重复出现时,即可省略繁琐的求和符号。

比如求和公式:

其中变量 a 和变量 b 的下标重复出现,则可将其表示为:

同理

进一步地,我们可以得到矩阵乘法的一个抽象

2. einsum的用途

einsum方法正是利用了爱因斯坦求和简洁高效的表示方法,从而可以驾驭任何复杂的矩阵计算操作。基本的框架如下:

以计算下式为例

可以写成如下形式

C = einsum('ij,jk->ik', A, B)

上述操作表示矩阵A与矩阵B的点积。输入的参数分为两部分:

  • 前面表示计算操作的指令串,
  • 后面是以逗号隔开的操作对象(数量需与前面对应)。

其中在计算操作表示中,

  • "->" 左边是以逗号隔开的下标索引,重复出现的索引即是需要爱因斯坦求和的;
  • "->" 右边是最后输出的结果形式。

在矩阵之间的运算中,下标可以分为两类:

  • 自由标(Free index),也就是在输入和输出端都出现的下标
  • 哑标(Summation index),在输入端出现但输出端没有出现的下标

实际上,即使是这个简单的例子,我测试使用 einsum 计算的速度要是常规手段的 3 倍。

用图像表示上述过程

A = np.array([[1, 1, 1],[2, 2, 2],[5, 5, 5]])B = np.array([[0, 1, 0],[1, 1, 0],[1, 1, 1]])

我们将轴标签的输入/输出轴画出来:

要理解上述计算过程,记住下面这三条规则

  • 输入arrays沿着重复的那个轴做乘法计算。

在本例中,轴标签的输入参数ij,jk中 'j' 被重复用了 2 次:1次表示A的第二个轴,1次表示B的第一个轴。这意味着要计算A的第二个轴(行方向)与B的第一个轴(列方向)的乘积。所以,此时要保证输入合法,就需要保证A的行长与B的列长一致。

  • 如果某个轴标签在输出标签中消失了,则表示要沿着该轴做求和计算。

此处,'j' 没有出现在输出标签中,这表示在执行乘法运算后,需要再沿着'j'轴做求和运算,求和运算减少了输出 array 的 1 个纬度。如果我们将输出标签改为'ijk',那么因为少了求和运算,最终得到的输出 array 将会是 3x3x3 形状的。

如果此处我们将输入/输出标签改为'ij,jk->',也即令所有的标签都不出现在输出轴标签里,那么我们将得到一个标量,这个标量是所有元素的和。

  • 我们可以获取任意顺序轴的结果

如果我们在轴标签中不写'->',那么 numpy 会将只出现一次的标签按照字母顺序组合,作为输出轴标签,所以 ij,jk->ikij,jk 效果上是等价的。指定轴顺序的输出,可以通过指定轴标签的顺序获得。例如,'ij,jk->ki'可以得到'ij,jk->ik'的转置矩阵。

了解上面三条基本规则后,再来看einsum如何计算矩阵乘法的就简单了。下图是左边是计算np.einsum('ij,jk->ijk', A, B)的结果,右图则是按照'j'轴求和后的结果:

按照前文所述,einsum是非常节约空间的计算函数,所以对于np.einsum('ij,jk->ik', A, B)einsum并不构建临时的 3D array 然后求和,而是直接在 2D 空间累加得到最终结果。

3. einsum的简单操作

下面两张表描述了 einsum 的基本操作。

若 A 和 B 为两个 1D arrays(假设在相应的操作中,A和B的形状总是合适的),那么:

Call signature NumPy equivalent Description
('i', A) A A
('i->', A) sum(A) A的所有元素和
('i,i->i', A, B) A * B A和B逐元素乘积
('i,i', A, B) inner(A, B) A和B的内积
('i,j->ij', A, B) outer(A, B) A和B的外积

若 A 和 B 为两个 2D arrays(假设在相应的操作中,A和B的形状总是合适的),那么:

Call signature NumPy equivalent Description
('ij', A) A A
('ji', A) A.T A的转置
('ii->i', A) diag(A) A 的对角
('ii', A) trace(A) A的迹
('ij->', A) sum(A) A的所有元素和
('ij->j', A) sum(A, axis=0) A的沿着axis=0的和
('ij->i', A) sum(A, axis=1) A的沿着axis=1的和
('ij,ij->ij', A, B) A * B A和B逐元素乘积
('ij,ji->ij', A, B) A * B.T A 和 B.T 逐元素乘积
('ij,jk', A, B) dot(A, B) A 和 B 的矩阵乘法
('ij,kj->ik', A, B) inner(A, B) A 和 B 的内积
('ij,kj->ikj', A, B) A[:, None] * B A的每一行与B的乘积
('ij,kl->ijkl', A, B) A[:, :, None, None] * B A的每一个元素与B的乘积

4. einsum轴标签中的'...'符号

在处理比较多的纬度时,为了方便,可以像 numpy array 一样使用 '...' 符号省略一些纬度的显式表示。例如,

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

这一行代码计算的是 a 的后两个轴与 2D array b 的乘积。

5. 注意事项

einsum 在求和时,不会提升数据类型。如果我们使用了位宽比较小的数据类型,可能会得到不期望的结果:

>>> a = np.ones(300, dtype=np.int8)
>>> np.sum(a) # correct result
300
>>> np.einsum('i->', a) # produces incorrect result
44

另外,einsum在 numpy 的计算库中并不一定总是最快的。类似于 dotinner 的函数一般链接到快速计算库 BLAS,可能会快于 einsum

参考文献

einsum方法详解(爱因斯坦求和) - 知乎

Python计算库NumPy的einsum(爱因斯坦和)使用简介 - 刘冲的博客

NumPy中einsum使用笔记相关推荐

  1. Numpy库中einsum函数用法

    Numpy中einsum函数用法 一.一维张量收缩 二.二维张量收缩 2.1 收缩到零维张量 2.2 收缩到一维张量 三.三维张量收缩(重难点) 3.1 例1 3.2 例2 四.其他功能介绍(次要) ...

  2. [云炬python学习笔记]Numpy中内置函数min(),max(),sum()与Python中内置函数min(),max(),sum()性能对比分析

    众所周知,Python有许多内置函数(例如min(),max(),sum()),Numpy也有自己的内置函数(np.min(),np.max(),np.sum()).由于Numpy的函数是在编译码中执 ...

  3. python科学计算笔记(一)NumPy中ndarray对象、ufunc运算、矩阵运算

    标准安装的Python中用列表(list)保存一组值,可以用来当作数组使用,不过由于列表的元素可以是任何对象,因此列表中所保存的是对象的指针.这样为了保存一个简单的[1,2,3],需要有3个指针和三个 ...

  4. python笔记之NUMPY中的掩码数组numpy.ma.mask

    python科学计算_numpy_线性代数/掩码数组/内存映射数组 1. 线性代数 numpy对于多维数组的运算在默认情况下并不使用矩阵运算,进行矩阵运算可以通过matrix对象或者矩阵函数来进行: ...

  5. numpy中matmul的使用(个人笔记)

    numpy中matmul的使用 简介:        numpy.matmul 函数返回两个数组的矩阵乘积.当两个数组都是二维数组的时候,就是数学上的两个矩阵的乘积. 例如: import numpy ...

  6. Python精讲Numpy基础,大牛笔记详细解释

    https://www.toutiao.com/a6664936105076326920/ 总认为Numpy是渣渣,直到深入接触以后才知道功能这么强大.堪比Matlab啊.果然是人生苦短,我用Pyth ...

  7. matlab imcrop 对应python函数_Python精讲Numpy基础,大牛笔记详细解释

    总认为Numpy是渣渣,直到深入接触以后才知道功能这么强大.堪比Matlab啊.果然是人生苦短,我用Python.所以本文作为一个记录&笔记,文章内容大多数取自网络以&官网快速入门等, ...

  8. python中把输出结果写到一个文件中_Python3.6笔记之将程序运行结果输出到文件的方法...

    Python3.6笔记之将程序运行结果输出到文件的方法 更新时间:2018年04月22日 14:27:32 投稿:jingxian 下面小编就为大家分享一篇Python3.6笔记之将程序运行结果输出到 ...

  9. pythonisnan_python - 在NumPy中快速检查NaN

    python - 在NumPy中快速检查NaN 我正在寻找最快的方法来检查NumPy数组np.nan != np.nan中NaN(np.nan in X)的出现.np.isnan(X)是不可能的,因为 ...

最新文章

  1. 博客摘录:网络管理员的两天
  2. 初学Java ssh之Spring 第三篇
  3. 通过几个问题深入分析Vue中的diff原理
  4. centos7配置静态IP
  5. Android中五种常用对话框的使用
  6. 打开黑色_表哥出差带回来一箱苹果,打开后发现是黑色的,大家表示都没见过...
  7. 如何判断文本文件的编码格式?
  8. spring.profiles.active配置了没生效_微服务架构之「 配置中心 」
  9. PS2018学习笔记(30-35节)
  10. insert sort java_java插入排序 Insert sort实例
  11. mount远程驱动器
  12. 通俗易懂的极限学习机(Extreme Learning Machine)
  13. 计算机考试祝福,考试前说的祝福语汇编35句 参加考试前的祝福语
  14. 一款免费好用的在线高效作图工具
  15. Word里面文字怎么加边框
  16. mysql导出到excel方法汇总
  17. Kienct与Arduino学习笔记(1) 基础知识之Arduino’Kinect‘Processing
  18. 三款截图软件:Snipaste+FastStone-Capture+FireShot
  19. Mybatis(一)——【快速入门、增删查改操作、核心配置文件描述及API】
  20. [error] Vivado代码仿真时错误提示:ERROR: [Common 17-39] ‘launch_simulation‘ failed due to earlier errors.

热门文章

  1. 从零开始在windows下使用QT根据点绘制图像
  2. bp神经网络的训练过程,BP神经网络图像识别
  3. Windows更新 (KB3045318)过程会暂停sql server相关服务
  4. 【负荷预测】长短期负荷预测(Matlab代码实现)
  5. WeCode在线少儿编程|创交会三大领域机器人各显神通
  6. 有限元求解两点边值问题之一
  7. 西门子SINUMERIK数控系统操作维修教学装置,QY-SKC25
  8. 如何布署 PyTorch 深度学习模型
  9. MR的原理和运行流程
  10. 最新深度技术GHOST XP系统旗舰增强版 V2016年