[转]矩阵运算优化问题

form https://www.zhihu.com/question/19706331

(一) MKL库问题

作者：过拟合
链接：https://www.zhihu.com/question/19706331/answer/25255444
来源：知乎
著作权归作者所有，转载请联系作者获得授权。

MATLAB的矩阵计算使用的是Intel自己出的Math kernel library（MKL）,这个库远比其他的blas/lapack库要快。C快在循环，要想矩阵计算也和MATLAB一样快，那就得链接MKL，写起来免不了各种折腾。而且，即使你链接上了，编译时各种优化选项之类的还是比不上人家专业的设定，速度很难接近MATLAB。
我自己在Gentoo上试过源里的所有blas/lapack库，无一能与MKL匹敌，而且连接近都不可能。甚至我把python的NumPy库链接上MKL后，速度也只是勉强接近。由于Gentoo的MKL库永远是最新的，而每一个新版本的MKL库对矩阵计算都有略微提升，导致可能暂时NumPy与MATLAB可以匹敌。但是一旦更新版本的MATLAB出来后，它会使用上更新的MKL库，这种领先优势就又丧失殆尽。你可以在MATLAB文档搜索中输入MKL，这样会被定位到MATLAB release notes，而里面就会含有这么一句话“Upgrade to Intel Math Kernel Libraries”，这就是每一个版本MATLAB矩阵计算都越发变态快的原因。
当然，刚我提到的python，其矩阵计算速度虽然微微落后于MATLAB，但是在很多其他地方是可以大大强于MATLAB的。例如绘制大规模三维点云，以及轻松调用gpu之类的。因此python在矩阵计算的微小速度劣势完全可以忽略，可以考虑用于科学计算。关于NumPy链接MKL，我之前写过一篇博文：Numpy使用MKL库提升计算性能

 

（二）MATLAB的优化

作者：王洋子豪
链接：https://www.zhihu.com/question/19706331/answer/16134037
来源：知乎
著作权归作者所有，转载请联系作者获得授权。

矩阵乘法是一个相对成熟的问题，根据矩阵的稀疏程度有不同的优化算法。
不使用GPU加速的MATLAB版本采用的是BLAS中的General Matrix Multiplication[1]。学术界有各种矩阵乘法算法将其复杂度降低到O(n^2.x)，例如Strassen和Winograd算法，在BLAS中应该已经使用了Strassen算法。
如果你的MATLAB是安装了Parallel Computing Toolbox的话，那么很可能它已经在使用GPU进行计算了。这种情况下采用的是MAGMA[2]。我没有使用过MAGMA，但我猜测它应该使用了cuBLAS来计算矩阵乘法。
宏观角度上对矩阵乘法的优化包括对局部内存使用的优化(Blocked/Tiled)以及对中间运算步骤的优化(Strassen/Winograd)，实现细节上的优化就非常繁多了。比如loop unrolling，多级的tiling，指令级并行等等。其中会牵扯到一些编译器和体系结构的知识，似乎对仅仅希望使用矩阵乘法函数的用户来讲没有什么太大必要去探究。所以我认为写出比MATLAB更快的矩阵乘法是可行但困难的。
[1] General Matrix Multiply
[2] MAGMA

 

（三）矩阵乘法cache优化

作者：杨树下的狐狸
链接：https://www.zhihu.com/question/19706331/answer/15392676
来源：知乎
著作权归作者所有，转载请联系作者获得授权。

偏个题。
说到这个想起之前朋友和我说到他最近在上一个课。
那个课上，教授要求他们写Cache friendly code。
尤其像矩阵这种很大块的东西，在运算时，会导致cache根本无法完全装下需要使用的数据。
此刻，如果程序没有设计得很有技巧，不断地刷新cache，会需要浪费大量的时间。
所以，他们教了好几种方法去计算矩阵，让整个计算过程中尽量减少cache的重新载入：
以下是引用朋友给我的邮件，作者是@Tian Tan :

先给你说个好玩的。这是我上的一门课的内容，叫high performance computing

这周在超级计算机上做了一个实验，
实验内容是想尽办法优化很小一段代码，比如矩阵乘法。

先说说cache的特点。
当访问内存中一个element的时候，
cpu会把那个element放到cache里面，
同时还会把它临近的elements放进cache。

衡量CPU快慢的一个标准是MFLOPS, 全称为millions of floating point operations per second. 实验的宗旨是写cache friendly code. 直接上例子吧。

A, B, C都是浮点数矩阵
for (i=0; i<N; i++)
for (j=0; j<N; j++)
for (k=0; k<N; k++)
C[i][j] += A[i][k]*B[k][j];

这是个很简单的矩阵乘法算法。但是这么写效率是不高的，
原因在于当N很大的时候，比如2048， 4096，会产生cache conflict。要解释这个术语比较麻烦，
想知道的话去看computer system: a programmer's perspective那本书的第六章。

但是矩阵乘法有别的写法。比如

for (i=0; i<N; i++)
for (k=0; k<N; k++)
for (j=0; j<N; j++)
C[i][j] += A[i][k]*B[k][j];

这种写法交换了k和j的位置，效率应该会比前面那个高些。（
术语是loop permutation）

还有的写法叫loop tiling，
tiling的实质是将大矩阵的乘法变成小的分块矩阵的乘法。

就用上一个例子吧。
for(it=0; it<N; it+=T)
for(jt=0; jt<N; jt+=T)
for(kt=0; kt<N; kt+=T)
for(i=it; i<it+T; i++)
for(j=jt; j<jt+T; j++)
for(k=kt; k<kt+T; k++)
C[i][j] += A[i][k]*B[k][j];

其中的T叫tiling size，能整除N。

这样先算小矩阵的话，cache 就能装下参与运算的elements, 对速度提升很大。

在实验中有一道题，
经过优化之后把运行时间从49秒降到13秒了。

矩阵乘法只是最简单例子，不同的code优化方式各异，
但是基本思想一样。

另外，正如上面说的。
使用针对你自己的CPU的编译器，编译器有可能能够识别到你的功能，做出相应的优化。