c++ - 如何在 C++ 中针对 x86 优化这三个矩阵的乘积？

Question

我有一个关键算法，其中大部分运行时间都花在计算密集矩阵乘积上:

A*A'*Y, where: A is an m-by-n matrix, 
               A' is its conjugate transpose,
               Y is an m-by-k matrix

Typical characteristics:
    - k is much smaller than both m or n (k is typically < 10)
    - m in the range [500, 2000]
    - n in the range [100, 1000]

基于这些维度，根据matrix chain multiplication的教训问题，很明显，将计算结构化为 A*(A'*Y) 从操作数的角度来看是最优的。我当前的实现就是这样做的，并且仅通过强制表达式的关联性来提高性能是显而易见的。

我的应用程序是用 C++ 为 x86_64 平台编写的。我正在使用 Eigen线性代数库，带有 Intel's Math Kernel Library作为后端。 Eigen 能够使用 IMKL 的 BLAS 接口(interface)执行乘法，并且在我的 Sandy Bridge 机器上从 Eigen 的原生 SSE2 实现迁移到英特尔优化的、基于 AVX 的实现的提升也很重要。

但是，表达式 A * (A.adjoint() * Y)(用 Eigen 的说法编写)被分解为两个通用矩阵矩阵乘积(调用 xGEMM BLAS routine)，中间创建了一个临时矩阵。我想知道，通过一次评估整个表达式的专门实现，我是否可以获得比我现在拥有的通用实现更快的实现。一些让我相信这一点的观察是:

• 使用上述典型维度，输入矩阵 A 通常无法放入缓存。因此，用于计算三矩阵乘积的特定内存访问模式将是关键。显然，避免为部分积创建临时矩阵也是有利的。

• A 及其共轭转置显然具有非常相关的结构，可以利用该结构改进整个表达式的内存访问模式。

是否有任何标准技术可以以缓存友好的方式实现这种表达式？我发现的大多数矩阵乘法优化技术都是针对标准 A*B 情况，而不是更大的表达式。我对问题的微优化方面很满意，例如转换为适当的 SIMD 指令集，但我正在寻找任何引用资料，以尽可能以对内存最友好的方式分解此结构。

编辑:根据目前收到的回复，我觉得我上面的内容有点不清楚。从我对这个问题的角度来看，我使用 C++/Eigen 的事实实际上只是一个实现细节。 Eigen 在实现表达式模板方面做得很好，但不支持将此类问题评估为简单表达式(仅支持 2 个一般稠密矩阵的乘积)。

在比编译器如何评估表达式更高的层次上，我正在寻找复合乘法运算的更有效的数学分解，旨在避免不必要的冗余内存访问由于 A 的共同结构及其共轭转置。结果可能难以在纯 Eigen 中有效实现，因此我可能只是在具有 SIMD 内在函数的专门例程中实现它。

Answer 1

这不是一个完整的答案(目前 - 我不确定它会成为一个完整的答案)。

让我们先稍微考虑一下数学。由于矩阵乘法是关联的，我们可以做(A*A')Y 或 A(A'*Y)。

(A*A')*Y 的浮点运算

2*m*n*m + 2*m*m*k //the twos come from addition and multiplication

A*(A'*Y) 的浮点运算

2*m*n*k + 2*m*n*k = 4*m*n*k

因为 k 比 m 和 n 小得多，所以很清楚为什么第二种情况要快得多。

但是通过对称性，原则上我们可以将 A*A' 的计算次数减少两次(尽管这对于 SIMD 可能不容易做到)因此我们可以减少 (A*A' )*Y 到

m*n*m + 2*m*m*k.

我们知道 m 和 n 都大于 k。让我们为 m 和 n 选择一个名为 z 的新变量，找出情况一和情况二相等的地方:

z*z*z + 2*z*z*k = 4*z*z*k  //now simplify
z = 2*k.

只要 m 和 n 都大于 k 的两倍，第二种情况的浮点运算就会更少。在您的情况下，m 和 n 均大于 100，而 k 小于 10，因此情况二使用的浮点运算要少得多。

在高效代码方面。如果代码针对缓存的有效使用进行了优化(如 MKL 和 Eigen)，那么大型密集矩阵乘法是计算限制而不是内存限制，因此您不必担心缓存。 MKL 比 Eigen 快，因为 MKL 使用 AVX(现在可能还有 fma3？)。

我不认为你能够比你已经使用第二种情况和 MKL(通过 Eigen)更有效地做到这一点。启用 OpenMP 以获得最大的 FLOPS。

您应该通过将 FLOPS 与处理器的峰值 FLOPS 进行比较来计算效率。假设您有一个 Sandy Bridge/Ivy Bridge 处理器。峰值 SP FLOPS 是

frequency * number of physical cores * 8 (8-wide AVX SP) * 2 (addition + multiplication)

对于双进动除以二。如果您有 Haswell，而 MKL 使用 FMA，则将峰值 FLOPS 加倍。要获得正确的频率，您必须对所有内核使用涡轮增压值(它低于单个内核)。如果您没有对系统进行超频，或者在 Windows 上使用 CPU-Z 或在 Linux 上使用 Powertop(如果您有超频系统)，则可以查看此内容。