线性代数计算是计算机的一个重要应用领域,进行线性代数计算的程序库一般使用标准的BLAS接口。其中矩阵乘法计算在BLAS中叫做GEMM(general matrix multiplication)。计算公式为:

C = alpha*A*B + beta*C

其中A为M×K矩阵,B为K×N矩阵,C为M×N矩阵,alpha和beta为系数。

GEMM也是机器学习领域的一个核心运算。深度学习模型中的全连接层和卷积层都可以用GEMM实现。GEMM运算量很大,与M×K×N成正比,所以高效的GEMM实现是很重要的,一般都是硬件厂商专门养活一帮人来写针对硬件优化的库,如NVIDIA的cuBLAS。然后写框架的来调这些库,然后写包的来调框架,然后我们调包。

这里研究一下如何在CPU上进行GEMM计算。

现代多级缓存CPU架构上的GEMM计算算法,可以参考Goto等人的论文Anatomy of High-performance Matrix Multiplication。一个思想就是分块,把矩阵(M, Matrix)分为条(P, Panel),再把条分为块(B, block),将GEMM转化为GEBP。

为什么要分块呢?是为了充分利用缓存,重用缓存中的数据进行计算。如果逐个计算的话,计算C矩阵的每一个元素都要读取A中的一行和B中的一列,计算下一个元素就要读取新的一行或者另一列,把原来放入缓存中的数据扔掉了,而读写内存操作的时间比计算操作要慢何止百倍。分块后就先把矩阵的一块读进来,用这些数据做更多的计算,提高计算操作与内存读写操作的比例。

GEBP的计算依然需要分块。一个GEBP计算,C=A×B,A可以一直留在缓存中,B和C一块一块的进入缓存进行计算。 由于现代CPU有多级缓存体系,从寄存器到L1/L2/L3缓存到RAM,所以具体实现更为复杂,需要考虑分块放在哪一级的缓存。而除了缓存层级的分块(MC, KC, NC),还有寄存器层级的分块(MR, NR),这是最小的计算单元了。

本文简单的介绍了一下矩阵乘的计算方法,理论上的具体细节可以直接看上面提到的论文。接下来看一看实际的工程代码中矩阵乘是怎么实现的。