8

我想找到一种尽可能快的方法将两个小布尔矩阵相乘，其中小表示 8x8、9x9 ... 16x16。这个套路会用到很多，所以它需要非常高效，所以请不要建议直截了当的解决方案应该足够快。

对于 8x8 和 16x16 的特殊情况，我已经有了相当有效的实现，基于此处找到的解决方案，我们将整个矩阵分别视为一个uint64_t或uint64_t[4]。在我的机器上，这比直接实现快大约 70-80 倍。

但是，在 8 < k < 16 的情况下，我真的不知道如何利用任何合理的表示来实现上述巧妙的技巧。

所以基本上，我愿意接受任何使用任何类型的表示（矩阵的）和函数签名的建议。您可以假设这针对 32 位或 64 位架构（选择最适合您的建议的）

4

4 回答 4

7

给定两个 4x4 矩阵 a= 0010,0100,1111,0001, b=1100,0001,0100,0100，首先可以计算转置 b' = 1000,1011,0000,0100。

然后得到的矩阵 M(i,j)=axb mod 2 == popcount(a[i]&b[j]) & 1; // 或奇偶校验

从中可以注意到，只要位向量适合计算机字，复杂度只会增加 n^2。

如果有一些特殊的置换和位选择操作可用，这至少可以加快 8x8 矩阵的速度。一个向量中的 NxN 位可以精确地迭代 N 次。（所以 16x16 几乎是极限）。

每个步骤由累加组成，即 Result(n+1) = Result(n) XOR A(n) .& B(n)，其中 Result(0) = 0，A(n) 是 A <<< n，并且 ' <<<' == 元素的列旋转，其中 B(n) 从矩阵 B 复制对角线元素：

    a b c          a e i          d h c          g b f
B=  d e f  B(0) =  a e i  B(1) =  d h c   B(2) = g b f
    g h i          a e i          d h c          g b f

在进一步考虑之后，更好的选择是^^^（逐行旋转）矩阵 B 并选择 A(n) == 列从 A 复制的对角线：

    a b c         a a a           b b b           c c c 
A=  d e f  A(0) = e e e , A(1) =  f f f,  A(2) =  d d d 
    g h i         i i i           g g g           h h h

编辑为了使以后的读者受益，我提出了便携式 C 中 W<=16 位矩阵乘法的完整解决方案。

#include <stdint.h>
void matrix_mul_gf2(uint16_t *a, uint16_t *b, uint16_t *c)
{
    // these arrays can be read in two successive xmm registers or in a single ymm
    uint16_t D[16];      // Temporary
    uint16_t C[16]={0};  // result
    uint16_t B[16];  
    uint16_t A[16];
    int i,j;
    uint16_t top_row;
    // Preprocess B (while reading from input) 
    // -- "un-tilt" the diagonal to bit position 0x8000
    for (i=0;i<W;i++) B[i]=(b[i]<<i) | (b[i]>>(W-i));
    for (i=0;i<W;i++) A[i]=a[i];  // Just read in matrix 'a'
    // Loop W times
    // Can be parallelized 4x with MMX, 8x with XMM and 16x with YMM instructions
    for (j=0;j<W;j++) {
        for (i=0;i<W;i++) D[i]=((int16_t)B[i])>>15;  // copy sign bit to rows
        for (i=0;i<W;i++) B[i]<<=1;                  // Prepare B for next round
        for (i=0;i<W;i++) C[i]^= A[i]&D[i];          // Add the partial product

        top_row=A[0];
        for (i=0;i<W-1;i++) A[i]=A[i+1];
        A[W-1]=top_row;
    }
    for (i=0;i<W;i++) c[i]=C[i];      // return result
}
于 2013-01-29T12:24:00.340 回答
5

如何将其填充到下一个“聪明”（例如 8 或 16）大小，对角线上的所有“1”？

于 2013-01-29T10:45:53.797 回答
4

根据您的应用程序，将矩阵及其转置存储在一起可能会有所帮助。您将节省大量时间，否则这些时间将在矩阵乘法期间用于转置，但代价是一些内存和更多操作。

于 2013-01-29T12:41:30.833 回答
1

有一种更快的方法可以使用 64 位乘法以及一些简单的位技巧来将 8x8 矩阵相乘，该方法适用于 GF[2] 或布尔代数。假设这三个矩阵在一个 64 位 int 中被打包成 8 位连续的 8 行，我们可以使用乘法来分散这些位并在一个 for 循环中完成这项工作：

uint64_t mul8x8 (uint64_t A, uint64_t B) {

    const uint64_t ROW = 0x00000000000000FF;
    const uint64_t COL = 0x0101010101010101;

    uint64_t C = 0;

    for (int i=0; i<8; ++i) {
        uint64_t p = COL & (A>>i);
        uint64_t r = ROW & (B>>i*8);
        C |= (p*r); // use ^ for GF(2) instead
    }
    return C;
}

如果您负担得起阻塞行以提高效率，则 16x16 的代码很简单。这个技巧也广泛用于高性能线性代数库，包括将矩阵划分为 N/M x N/M 个 MxM 子矩阵块，选择 M = 2^m 以最大化缓存中的局部性。处理 N % M != 0 的常用方法是用 0 填充行和列，以便可以对所有块乘法使用相同的算法。

我们可以将相同的想法应用于可变维数 8 >= N >= 16 的布尔矩阵，只要我们能够以行阻塞格式在内部表示矩阵。我们只是假设矩阵是 16x16，最后 16-N 行和列用 0 填充：

void mul16x16 (uint64_t C[2][2], const uint64_t A[2][2], const uint64_t B[2][2]) {

    for (int i=0; i<2; ++i)
        for (int j=0; j<2; ++j)
            C[i][j] = mul8x8(A[i][0],B[0][j])
                    | mul8x8(A[i][1],B[1][j]); // once again, use ^ instead for GF(2)
}

请注意，我们仅在 8x8=64 整数乘积和一些位运算中完成了 16x16 矩阵乘法。

mul8x8 也可以通过现代 SSE/AVX 矢量指令得到很大改进。理论上可以用一条 AVX512 指令并行执行所有 8 个乘积（我们仍然需要先将数据分散到 ZMM 寄存器），然后使用 lg2(8) = O(3) 指令进行水平归约。

于 2019-03-22T20:45:54.587 回答