关于高性能数值运算

转自林达华的博客 dahuasky.wordpress.com （需翻墙）

经过一个月的繁忙后，在Microsoft Research的工作进入收尾阶段。于是，也有时间上来写点东西了。在微软intern期间有很多的收获，等回到学校后，回顾整个过程，再写个总结吧。这篇blog里，主要说的是一个实用的东西，就是在一般的程序设计中进行高性能数值运算的问题。

数值计算是我们实现算法过程中必不可少的部分，在大规模计算中（比如在大型数据集上训练模型，或者在多节点的网络中进行统计推断）对计算性能的追求是非常必要的（在关键环节的优化甚至能把实验时间缩短10倍以上）。对这个问题的思考其实事出有因。一直以来，在学校里面都是使用matlab作为主要的计算工具，可是，到了微软之后情况发生了很大的变化——

在外面看来，微软虽然是很有钱的公司，不过并不是对什么事情都这么慷慨的（个人觉得MIT在不少地方比微软慷慨得多～～嗯，除了工资）。在Redmond的研究院有几百名正式的研究员，外加几百intern；而matlab的license只有30个，换言之，就是整个研究院范围只能同时开30个matlab进程。每天在公司的邮箱里，总能收到很多请求别人release matlab license的邮件，语气之恳切近乎哀求。在这样的条件下，我的mentor跟我说：matlab is simply not an option。在这种情况下，必须在脱离matlab的环境中寻求高性能计算的其它途径。

从MATLAB矩阵乘法到BLAS dgemm

事实上，无论是在matlab，还是在别的语言下的计算库，高性能计算过程都大概包含了三个层次。为了说明问题，还是以matlab中的矩阵乘法为例子说起。

在matlab中，我们要计算两个矩阵 A 和 B 的乘积，我们会在命令行或者m-file里面写 A * B。这么一个简单的语句是如何转化为实际的计算过程的呢？ 首先，matlab的解释器会分析这个语句的语法，从而知道你要做的是计算 A 和 B 的乘积。但是，matlab 其实只是个外包装，并不是自己进行这个计算的，而是把计算过程委托给一个核心的C语言计算函数叫做dgemm，然后dgemm进行实际的计算，最后matlab对dgemm输出的结果重新打包返回。

这个名字很奇怪的函数究竟是怎么回事呢？它其实是BLAS Level 3中最重要的函数 (常用BLAS/LAPACK 或者 MKL的朋友对这个函数肯定不陌生)。BLAS 全称 Basic Linear Algebra Subprograms，它是在1979年被提出的来用于支持建立一个叫做LAPACK的矩阵代数的运算库。到了今天，BLAS函数已经成为矩阵和向量运算的事实上的“国际标准”接口。BLAS的函数分为三个Level,

Level – 1是向量和向量之间的运算，比如点积 （ddot)， 加法和数乘 (daxpy), 绝对值的和 （dasum)， 等等；

Level – 2 是向量和矩阵的乘法运算，最重要的函数是一般的矩阵向量乘法(dgemv)

Level – 3 是矩阵和矩阵的乘法运算，最重要的函数是一般的矩阵乘法 (dgemm)

dgemm这个名字是什么意思呢？其实全称可以翻译为 double-precision generic matrix-matrix muliplication. 其实，Level-3虽然只做一件事情，就是矩阵相乘，但是，它还有很多其它的函数：比如sgemm (单精度一般矩阵乘法)，dsymm (双精度对称矩阵乘法)，zhemm(双精度复数埃米特矩阵乘法)，诸如此类还有很多 。。。。。之所以要分成这么多种，主要是针对每种不同类型的矩阵，都要分别针对性地设计专门的算法使得对它的性能尽可能的高。

另外还有一套同样著名的标准运算接口，叫做LAPACK(Linear Algebra PACKage)，它是基于BLAS的基础上建立的线性代数运算库，提供一些更复杂的功能，比如LU, QR, SVD分解，或者求解特征值和特征向量，又或者求解线性方程组以及最小二乘法等问题。

其实 BLAS 可以理解为只是一套函数接口，在遵循接口的情况下，不同的库可以有不同的实现。很多著名的软硬件厂商都分别实现了针对自己产品专门优化过的BLAS实现。最著名的实现有下面一些：

ATLAS: 这是一套开源的实现。很多开源的数值软件都使用ATLAS优化的BLAS；

Intel MKL：这是Intel针对它的CPU系列开发的核心数学运算库，提供了完整的BLAS和LAPACK实现，除此以外，还有很多功能；

AMD ACML： 作为Intel的长期对手，自然也不甘人后；

Sun Performance Library: 主要针对Sun的SPARC架构；

另外，Apple, HP, NEC都有各自的BLAS实现。各大芯片厂商其实都不遗余力地提高自己的BLAS实现的性能——因为，前面提到的那个著名函数dgemm的运算速度是衡量一个CPU数值运算性能的非常重要的指标。一般来说，我接到一台新的机器，都会打开matlab算几次1000*1000或者2000*2000的随机矩阵乘法（前面提到，matlab也是调用核心的dgemm做这个事情），这样我就大体知道这台机子的数值性能了。

既然大家都实现了BLAS，为什么我们通常不直接用它呢？答案很简单，太不方便了。dgemm这个函数做矩阵乘法，一般可能觉得它会有三个参数，两个输入，一个输出。考虑到需要把大小信息m, n, l 也提供进去，最多6个参数。好，我们一起来参观一下这个C函数的声明：

void cblas_dgemm(const enum CBLAS_ORDER Order,

const enum CBLAS_TRANSPOSE TransA,

const enum CBLAS_TRANSPOSE TransB,

const int M,

const int N,

const int K,

const double alpha,

const double *A,

const int lda,

const double *B,

const int ldb,

const double beta,

double *C,

const int ldc);

总共14个参数。不知道大家是否有心情在一般的算法开发中用这东西做矩阵乘法，反正我是不会这么做的。而matlab这类软件的好处是包装了这些很难用的函数，使得我们在很愉快的写A * B这样的表达式时，同时享用BLAS带来的高性能。

让矩阵乘法变得更快

我们都知道矩阵乘法的运算规则，其实很简单——一个10行内的的三重循环就能实现。为什么那么多大公司还要年复一年的研究怎么算矩阵乘法（还有别的更简单的数值运算过程）呢？因为，算对很容易，算得快却非常艰难。通常来说，会使用这么一些技术：

－ 分解小矩阵块，每个小块都用特殊优化过的算法计算；

－ 把循环解开，减省内重循环中用于更新循环变量的开销。对于矩阵乘法这样的高密度数值运算，内重循环中每次循环哪怕减少一个指令周期，对整体性能都是很大的提高；

－ 对高速缓存（Cache）的调度进行改进，提高命中率；

－ 充分利用指令流水线，使得数据的读取，计算，写入这些接续操作在流水线内同时进行；

－ 使SIMD指令集（通常是SSE, SSE2, SSE3等）。SIMD全称是Single Instruction Multiple Data，就是一条指令在单周期内处理多组数据，在新的CPU中都支持这些指令。

SIMD对于向量运算的加速有重要意义。举个例子，在SSE2有指令能在单周期同时对一组128位浮点数（两个double，或者四个float）进行加减乘除或者开平方运算，可以想见，充分利用这些指令能对计算速度成倍提高。对于不同的CPU，不同的架构，不同大小的矩阵，对这些技术的运用也不尽相同，因此最终形成的实现算法非常非常复杂——资料显示，现代CPU上高度优化的矩阵乘法实现不是10行以内的for-loop，而是超过10万行！其中包含了大量和CPU和Cache密切相关的东西。

归纳起来，高速矩阵运算的全过程包含三个层次：

1. 上层封装好的接口： 比如matlab里面的A * B;

2. 中间层的 BLAS/LAPACK 的 C 函数，像dgemm，它们通过非常复杂并且硬件相关的技术来提高速度。

3. 这些C函数会调用CPU中的指令进行运算，其中SIMD（SSEx)指令集对于运算速度的提高有着关键意义。

提高基本运算的速度显然不是我们的任务，Intel和AMD的工程师会管这个事情。对于我们来说，合理的做法就是拿来主意，让他们开发出来的技术来提高我们的实验效率。如果能用matlab，matlab已经把这个包装好了。如果你是那种对性能特别有追求的人，那么你可以到Intel或者AMD的网站获取最新版的MKL或者ACML，安裝到你的机器上，然后设置BLAS_VERSION环境变量就能让matlab使用最新的高性能库了。

其它的常用数值软件，比如Mathematica, numpy, Octave, R都提供了方法链接不同的BLAS实现。对性能关键的数值运算程序，在选用支持库时，能否接入高性能BLAS实现是一个重要的标准。如果两个库都使用
98ec
同一种BLAS实现，那么它们的计算性能应该是差不多的（除了在封装层的overhead略有差别）。而接口封装的质量当然也是重要的考虑，这直接关系到它的易用性——这方面matlab是我见过做得最好的，Octave的接口仿造matlab，不过octave的运行环境比matlab就差多了。numpy和R都存在类似问题。

－－－－－

至于我自己，在微软intern期间，按照自己的设计思路重新用C++设计了一套新的矩阵库（这不是我在微软的主要工作，只是磨刀不误砍柴工而已），主要的设计目标有两个：(1)向下为接入MKL之类的BLAS/LAPACK实现提供方便；(2)向上对用户 提供类似matlab的接口。同时也利用这个机会探索一些有趣的C++ template metaprogramming的技术。在微软写的代码自然是不能公开的。不过，我每天回到住处后又改进设计思路重写了一套（这套就不属于微软了），准备带回学校后继续使用。在经过一段时间考验后，我想可能会把源代码开放出来和大家交流。