Jas*_*n R 9 c++ optimization x86-64 matrix-multiplication eigen
我有一个密钥算法,其中大部分运行时用于计算密集矩阵产品:
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]
基于这些维度,根据矩阵链乘法问题的教训,很明显,在运算数意义上将计算结构化为最优A*(A'*Y).我当前的实现就是这样做的,而只是强迫关联性到表达式的性能提升是显而易见的.
我的应用程序是用C++编写的,用于x86_64平台.我正在使用Eigen线性代数库,英特尔的数学核心库作为后端.Eigen能够使用IMKL的BLAS接口来执行乘法,并且从我的Sandy Bridge机器上移动到Eigen的原生SSE2实现到Intel优化的基于AVX的实现的提升也很重要.
然而,表达式A * (A.adjoint() * Y)(用Eigen的说法)被分解为两个通用的矩阵 - 矩阵乘积(调用xGEMMBLAS例程),在它们之间创建一个临时矩阵.我想知道,通过一次专门的实现来一次评估整个表达式,我可以得到一个比我现在的通用更快的实现.一些让我相信这一点的观察结果是:
使用上述典型尺寸,输入矩阵A通常不适合缓存.因此,用于计算三矩阵乘积的特定存储器访问模式将是关键.显然,避免为部分产品创建临时矩阵也是有利的.
A 并且它的共轭转置显然具有非常相关的结构,可以利用它来改善整体表达的存储器访问模式.
是否有任何标准技术以缓存友好的方式实现这种表达式？我发现的矩阵乘法的大多数优化技术都是针对标准A*B情况而不是更大的表达式.我对这个问题的微优化方面很满意,例如转换成适当的SIMD指令集,但是我正在寻找任何可以用尽可能最友好的方式打破这个结构的引用.
编辑:根据到目前为止的反应,我想我上面有点不清楚.我使用C++/Eigen的事实实际上只是我对这个问题的看法的一个实现细节.Eigen在实现表达式模板方面做得很好,但是不支持将此类问题作为简单表达式进行评估(仅支持2个通用密集矩阵的产品).
在比编译器如何评估表达式更高的层次上,我正在寻找复合乘法运算的更有效的数学分解,并且由于共同的结构A及其共轭转置而倾向于避免不必要的冗余存储器访问..结果可能很难在纯Eigen中有效地实现,所以我可能只是在具有SIMD内在函数的专门例程中实现它.
这不是一个完整的答案（但我不确定它会成为一个完整的答案）。
让我们先考虑一下数学。由于矩阵乘法是关联的，我们可以执行 (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，因此情况 2 使用的浮点运算要少得多。
在高效的代码方面。如果代码针对缓存的有效使用进行了优化（如 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。