矩阵乘法的顺序：同一个算法快几十倍的秘密
前言
起因是同事说同样三重循环算矩阵乘法，只是把循环变量的顺序换了一下，快了几十倍，一开始我是不信的，算法都没换哪来的几十倍。自己跑了一遍被啪啪打脸，顺手把背后的原理捋清楚了，故有此文备忘。
先看现象
C语言里教科书式的矩阵乘法，N阶方阵要做 $2N^3$ 次浮点运算，一行注释都不用多加：
// C = A * B, N阶方阵，row-major存储
void gemm_ijk(double *C, const double *A, const double *B, int N) {
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
C[i*N+j] += A[i*N+k] * B[k*N+j];
}只把j和k这两层循环交换一下位置，语义完全等价：
void gemm_ikj(double *C, const double *A, const double *B, int N) {
for (int i = 0; i < N; i++)
for (int k = 0; k < N; k++)
for (int j = 0; j < N; j++)
C[i*N+j] += A[i*N+k] * B[k*N+j];
}在普通笔记本上拿N=2048实测，gcc开O2：
ijk: 约 9 秒
ikj: 约 1.5 秒一个字母顺序的事，六倍没了，规模再大差距还更明显。
为什么会这样
CPU算加减乘其实快的离谱，真正慢的是等数据从内存里搬过来。寄存器一拍就到，L1缓存几个周期，L2十几个周期，真掉到内存就是上百个周期。所以现代CPU的性能瓶颈往往不在计算而在访存，矩阵乘法刚好是个重灾区。
关键在C语言二维数组是按行连续存的（row-major），而CPU搬内存不是按需单个搬，是一次搬一整条cache line（通常64字节），搬过来的同时还会把相邻的数据一起请进缓存。
再看上面两个循环：
- ijk版：最内层沿着j走，每次访问的是B[k*N+j]，这个是沿行访问没问题。但回头看，B的同一行在被不同的i反复用到时，中间隔了整整一层的k循环，缓存早就被冲得干干净净，基本每次都要重新去内存拿。
- ikj版：最内层固定住A的一个元素A[i*N+k]，然后沿着j一口气扫完C的一整行和B的一整行，这两行数据进了缓存就被最内层循环吃干榨净，访存局部性直接拉满。
一句话总结：同样的计算量，一个访存模式是把缓存当摆设，一个是把缓存喂饱。著名论文"Anatomy of HPL"里管这种分析叫roofline模型，用算术强度 $I = \dfrac{\text{FLOPs}}{\text{Bytes}}$ 衡量每搬一个字节能换来多少计算：算力是屋顶，内存带宽是屋顶下面的坡，强度太低跑在坡上，内存带宽跟不上，算力再高也白搭。
还能再快吗
能，而且快得多。ikj版本在矩阵大到缓存装不下整行整列的时候又会退化，于是有了分块（blocking/tiling）：把大矩阵切成一个个能塞进缓存的小块，先算块内的小矩阵乘法再拼起来，保证数据进一次缓存就被用够本再出去。
再往深处走还有寄存器分块、SIMD向量指令、汇编级调度，这一套组合拳打包起来就是BLAS库。同样的N=2048：
朴素ijk: 约 9 秒
OpenBLAS: 约 0.05 秒两个数量级的差距，这就是为什么所有正经的数值计算都要求你别手写矩阵乘法，直接调库。
BTW
当前热门深度学习，训练时说到底绝大部分时间都在做矩阵乘法，GEMM（通用矩阵乘）就是那个真正吃算力的角色。GPU上的Tensor Core、各种加速卡本质上都是把GEMM做到极致的专用电路，所以理解了这一篇，也就理解了为什么大模型烧的是显卡而不是CPU。