向量矩阵乘法、浮点向量、二进制矩阵

Posted 2023-02-16

【中文标题】向量矩阵乘法、浮点向量、二进制矩阵

【英文标题】：Vector matrix multiplication, float vector, binary matrix

【发布时间】：2019-10-14 15:58:09

【问题描述】：

我想将一个大小为 N 的浮点向量与一个大小为 NxM 的矩阵相乘。

矩阵是一个二元矩阵（只包含零和1），相对稀疏：非零值的密度在1-5%之间。

目前，我将其形成为密集向量和稀疏浮点矩阵乘法。

但这只是矫枉过正，不是吗？

如果我将矩阵的列存储为位集，然后乘法只是使用位集作为索引向量，然后将其相加，会怎样。

我假设我可以将其形成为 SSE/AVX 中的矢量化操作，例如 load + and + sum 或 load + mask + sum

如果你能指出我这样做的正确内在函数，我将不胜感激，主要问题是处理解包 bitset 的最佳方法是什么？

【问题讨论】：

如果你的密度永远不会超过 5%，你可能最好只存储非零的索引并执行vgatherdps（单个收集比负载更昂贵+屏蔽，但它不应该贵 20 倍）——当然要确保你有一些独立的累加器寄存器，以避免延迟。如果你的密度模式中有一些结构（例如，某种块结构），你的矩阵也可能有更好的表示。我想在你有更高的密度之前，bitset 的想法不会更好（不过我没有对此进行基准测试）。

好点。实际上，我正在检查可能发生的各种场景和我将要处理的矩阵，密度更接近 5-15%，我会更新问题。我现在担心的是内存使用情况，对于行列索引，每个非零单元格需要两个整数，对于大的 N 和 M，可能会很重要，N 大约是 100K 到 500K，M 大约是 500k 到 1M

使用压缩行存储，您应该只需要一个整数（大到足以容纳最高列号）每个（非零）元素和一个整数（大到足以容纳元素的数量) 每行。在经典 CRS 中，您还需要存储实际的非零值，但由于在您的情况下它们都是 1.0f，因此您不需要存储它们。即，只要你有少于大约 1/32 个非零值（对于 32 位索引），就记忆而言，你会更好。如果您的非零值有些聚集，则某些混合方法可能会起作用（例如，将位掩码块及其坐标存储在矩阵中）。

您还打算编辑/澄清您的问题吗？如果您的矩阵（部分）具有一些规则结构，那将特别有趣——但这可能最好是一个新问题。

【参考方案1】：

所以结果向量的每个元素都是输入向量的掩码总和？这些掩码来自矩阵的列，因此它们不是连续的位。

使用连续位图的屏蔽总和对于 AVX512 来说是微不足道的（只需使用合并屏蔽加法或零屏蔽加载）。对于 SSE/AVX2，您将使用 is there an inverse instruction to the movemask instruction in intel avx2? + _mm256_and_ps。或者对掩码向量进行优化的一些变化，例如具有 32 位广播负载，然后将其转移到下一步。而不是为每个字节进行另一个未对齐的 dword 广播。

但是你的掩码位 not 是连续的，你有一个选择：

分别处理每个输出向量元素，最后是水平和。需要收集位并制作矢量掩码。可能很难，除了 M=32 的情况，其中位跨度已经与连续的 32 位浮点数对齐。

累积 4 或 8 个输出元素的向量

，使用 4 或 8 个掩码位的连续组。因此，您在外循环上进行矢量化，在输入向量的内循环中进行广播加载。

使用这个。

您实际上应该

使用多个向量和展开以隐藏 FP 添加延迟。

像__m256 v = _mm256_set1_ps(invec[i]) 这样的广播负载基本上是免费的（vbroadcastss 是纯负载uop，没有ALU shuffle uop）。即使在循环结束时，您也不需要对浮点数进行任何其他改组，只需纯垂直 SIMD：您只需 _mm256_storeu_ps 进入输出向量。

而且您使用的是连续的掩码位组，因此通常的反向移动掩码问答很有用。

  // untested, rough outline of what it might look like

  uint8_t matrix[rows * cols];  // bit matrix in chunks of 8 bits
  float invec[N], outvec[N];    // A normal function will just take pointer inputs.

  constexpr int unroll = 4;
  for(int outpos = 0 ; outpos < M-8*unroll+1 ; outpos += 8 * unroll) 
      __m256 sum0, sum1, sum2, sum3;  //optionally use an array of accumulators, sums[unroll];
      sum0 = sum1 = sum2 = sum3 = _mm256_setzero_ps();
            // optionally peel the first inner iteration to just load+mask without adding to 0.0
      for (int inpos = 0 ; in < N ; in++ )
          __m256 inv = _mm256_set1_ps(invec[inpos]);
          __m256 mask0 = inverse_movemask(matrix[outpos*stride + inpos + 0]);  // 8 bits -> 8 vector elements
          __m256 mask1 = inverse_movemask(matrix[outpos*stride + inpos + 1]);
          ...

          sum0 = _mm256_add_ps(sum0, _mm256_and_ps(inv, mask0) );  // add in[i] or 0.0 according to mask
          sum1 = _mm256_add_ps(sum1, _mm256_and_ps(inv, mask1) );
          ...
      
      __m256_storeu_ps(&outvec[outpos + 0*8], sum0);
      __m256_storeu_ps(&outvec[outpos + 1*8], sum1);
      __m256_storeu_ps(&outvec[outpos + 2*8], sum2);
      ...
  

  not-unrolled __m256 and/or __m128 cleanup for M % (8*unroll) != 0

  cleanup for M % 4 != 0 using __m128 broadcast loads 
    for the last 1..3 rows of masks
    maybe use a masked store (AVX2 vmaskmov) or pad your output vector

每个内部循环迭代以不同的方式掩盖一个浮点 8 * unroll，并累积到相应的 8 * unroll 不同的运行总数中。（跨越 unroll 向量，每个向量有 8 个浮点数。）

---

这对内存带宽也很有用

在 vec*mat 乘积中，您只读取每个位图位一次，但输入向量被有效使用M 次。在连续的位图行上循环提供了良好的局部性，不需要多次加载任何这些缓存行。

即使使用 AVX512 和每个时钟 2x _mm512_mask_add_ps，每个 FP 元素添加 1 位对于位图加载来说带宽并不多。

但是，您循环输入向量 M/(8*unroll) 次。每个和向量的掩码加法使用不同的掩码位但相同的广播输入float。由于矩阵元素比向量元素小 32 倍，这还不错。

每 4x 或 8x vaddps 指令加载一个浮点数是非常好的计算强度。尤其是在没有 AVX512 的情况下，位图 -> 向量蒙版会花费周期。

为了在缓存/内存带宽方面提供更多帮助，cache-blocking / loop-tiling for L2 cache size (256kiB) 可能有助于重复使用输入向量元素。但我不确定您是否可以有效地阻止输入和输出。与 mat*mat 产品不同，只有 O(n^2) 工作要做。重新读取输入并仅写入一个输出流可能很好，但您可以找到一个中间立场，将部分结果添加到输出向量的部分块中。但是不再在一个连续的流中读取位矩阵。只要你停在缓存行边界就可以了。

---

如果您的NxM 矩阵恰好有M = 32，那么它与float 的大小完全匹配，并且_mm256_loadu_si256 将获得一个向量，该向量在每个矩阵的低位中都有outvec[0] 的掩码位元素。以及高位中outvec[31] 的掩码位。您可以使用_mm256_blendv_ps 将它们应用于和的输入，然后左移 1 以将下一位向上移动到顶部位置。 （vblendvps 的替代方案是 psrad by 31 + andps：算术右移将最高位广播到所有位置。

但这可能并不比其他方式更好，即使对于这种特殊情况也是如此。您可以展开不同向量中的多个输出元素，以便多次重复使用浮点向量。

---

使用 AVX512F，您可以将矩阵行用作 a masked add like _mm512_mask_add_ps 的 __mmask16 值。sum = _mm512_mask_add_ps(sum, matrix[col*rowstride + row], sum, invec[i]); 如果 matrix 是 uint16_t 的数组。

或使用 AVX512BW，kmovq 64 位掩码到 k 寄存器和 kshift 向下，以匹配展开 4 个矢量累加器。不幸的是，kmov k, [mem] 在 Skylake-X 上是 2 微指令：负载 + 端口 5，而不仅仅是可以写入掩码 reg 的负载微指令。因此，使用kshift 加载 3x 解包是纯粹的胜利，而 4x kmovw k1, [mem] / kmovw k2, [mem+2] 等。如果没有每个端口 5 uop，就无法在k 寄存器的底部获取每个 16 位掩码数据一。因此，它与具有 2 个 FMA 单元的 SKX 内核上的 512 位 FMA/add/mul 吞吐量竞争，否则只是前端吞吐量成本。

【讨论】：

+1 这太完美了！从文件中读取矩阵时，我确实可以将这些位形成为连续的。我忘了说这是在优化问题中使用的，并且有一个大矩阵（N 是 100K 到 500K，M 大约是 500K 到 1M），它在运行过程中根本不会改变。这个运行的平台是高度异构的，都有 AVX2 向上，很多没有 AVX512。那么你认为我应该采用“使用连续位图的掩码总和......”段落中的解决方案吗？

@NULL：那么别忘了投票。 :P 你可以制作一个 AVX512 版本并进行运行时调度； 100k x 500k 足以摊销调度开销。但是，是的，我强烈建议在每次传递给向量时计算几个 outvec[i + 0..7] 或 0..15 的 SIMD 向量。这足以考虑循环平铺（也称为缓存阻塞），比如只做适合 L2 缓存的输入向量的一部分，然后让后面的块将它们的部分结果添加到输出向量中，而不是仅仅存储。

@NULL：添加了循环应该是什么样子的粗略轮廓。

两个小问题： 1- 在内部循环中 _mm256_set1_ps 广播一个浮点数，但我们不应该用 8 个浮点数填充它吗？因为否则那将是多余的，不是吗？ inverse_movemask 也在做一个 8 位（所以矩阵中 col 的 8 个元素），那么 and+add 需要 8 个浮点数？ 2- 这里的步幅只是 N 对吗？

@NULL：它是一个浮点数，每个向量使用 8 个掩码位。请记住，我们在 outer 循环上展开/矢量化，每次迭代执行更多 outvec[i] 元素。每次内循环迭代以不同的方式屏蔽一个浮点数8 * unroll，并累加到对应的8 * unroll不同的运行总数中。步幅是N/8 向上舍入到一个完整字节，因此每一行位图都是字节对齐的。或者四舍五入，这样如果你愿意，每一行都是双字对齐的。