Eigen中矩阵-矩阵乘法比多次矩阵-向量乘法更慢的原因及优化

为什么多次矩阵-向量乘法反而更快？ #

这种现象主要源于Eigen在不同计算场景下的优化策略、内存访问模式以及并行开销的差异，具体原因包括：

• 内存访问与缓存命中率差异：
  Eigen默认采用列优先存储，矩阵-向量乘法（gemv）会连续访问矩阵的列与向量的元素，缓存局部性极佳，几乎没有缓存失效。而矩阵-矩阵乘法（gemm）虽然也利用缓存，但当Z的列数极少（比如3列）时，矩阵块的复用率很低，额外的内存寻址反而会降低缓存效率。

• 并行化的开销收益比失衡：
  Eigen的并行优化针对大规模矩阵乘法设计，启动线程、任务调度都有固定开销。当Z的列数很少时，并行计算带来的收益不足以抵消调度开销，反而不如单线程下的矩阵-向量乘法高效。而多次矩阵-向量乘法可以更灵活地利用CPU的指令级并行（ILP），循环展开、指令调度的优化更充分。

• 额外的数据复制开销：
  将u、v、w合并为矩阵Z，以及将结果D拆分回多个向量的过程，会产生额外的内存复制操作，这部分开销在计算量不大时会成为性能瓶颈。而直接执行矩阵-向量乘法，结果可以直接写入目标向量，无中间复制成本。

• Eigen的专用优化路径：
  Eigen对矩阵-向量乘法有专门的优化实现（如调用BLAS的gemv），针对向量计算的循环展开、寄存器利用更精细；而矩阵-矩阵乘法的gemm优化更偏向大矩阵的分块计算，小列数场景下无法发挥优势。

如何优化这类计算？ #

根据你的场景，可从以下几个方向优化：

• 保留多次矩阵-向量乘法的实现：
  如果最终需要的是多个独立向量结果，直接执行多次a_stack = M_stack * u这类操作即可，省去合并/拆分矩阵的额外开销，这也是当前最优的选择。

• 匹配存储顺序提升缓存效率：
  确保M_stack的存储顺序与计算模式匹配。例如，如果你的计算更偏向行访问，可显式指定行优先存储（Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>），减少缓存失效次数。

• 按需启用并行优化：
  若后续需要处理的向量数量大幅增加（比如数十列以上），可开启Eigen的并行支持：

  • 编译时添加-fopenmp并定义EIGEN_USE_OPENMP
  • 或使用MKL后端（定义EIGEN_USE_MKL_ALL）
    但注意向量数量较少时，关闭并行（定义EIGEN_DONT_PARALLELIZE）反而更快。

• 手动分块优化大矩阵：
  如果M_stack是超大规模矩阵（比如数万行/列），可手动将其拆分为固定大小的子块（如256x256），对每个子块执行矩阵-向量乘法，利用局部性原理最大化缓存利用率。

• 避免不必要的临时变量：
  直接在目标向量上进行计算，比如a_stack.noalias() = M_stack * u，使用noalias()避免Eigen创建临时矩阵，进一步减少内存开销。

内容的提问来源于stack exchange，提问作者Andrea Gotelli