如何在C++中为x86优化这三个矩阵的乘积

How to optimize this product of three matrices in C++ for x86?

本文关键字：三个 C++ 中为 x86 优化      更新时间：2023-10-16

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

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接口来执行乘法，从Eigen的原生SSE2实现到Intel在我的Sandy Bridge机器上优化的基于AVX的实现，这一提升也很重要。

然而，表达式A * (A.adjoint() * Y)(用特征语写成)被分解为两个通用矩阵矩阵乘积(对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，所以情况二使用的浮点运算要少得多。

就高效代码而言。如果代码是为了高效使用缓存而优化的(就像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，如果你有超频系统。

使用一个临时矩阵来计算a'*Y，但要确保告诉特征没有混叠：temp.noalias() = A.adjoint()*Y。然后计算结果，再次告诉本征对象没有别名：result.noalias() = A*temp。

只有在执行(A*A')*Y时才会有冗余计算，因为在这种情况下，(A*A')是对称的，只需要一半的计算。然而，正如您所注意到的，执行A*(A'*Y)仍然快得多，在这种情况下没有冗余计算。我确认，临时创建的成本在这里完全可以忽略不计。

我想执行以下

result = A * (A.adjoint() * Y)

将与相同

temp = A.adjoint() * Y
result = A * temp;

如果你的矩阵Y适合缓存，你可能会像一样利用它

result = A * (Y.adjoint() * A).adjoint()

或者，如果不允许以前的记法，比如

temp = Y.adjoint() * A
result = A * temp.adjoint();

然后你不需要做矩阵A的伴随，并存储A的临时伴随矩阵，这将比Y的矩阵昂贵得多。

如果你的矩阵Y适合缓存，那么在第一次乘法的a列上循环，然后在第二次多重乘法的a行上循环应该会快得多(第一次乘法在缓存中有Y.attaint()，第二次乘法在temp.account())，但我猜内部特征已经在处理这些问题了。

相关文章：