二进制矩阵乘法位操作技巧

Question

摘要

您好，假设您有两个不同的独立的64位二进制矩阵A和T（T是另一个以转置形式存储的矩阵，在乘法过程中使用矩阵的转置版本可以对T的行而不是列进行操作，这对于二进制算术非常方便），您想要将这些矩阵相乘。唯一需要注意的是矩阵乘积结果被截断为64位，如果在某些特定的矩阵单元格中得到大于1的值，则结果矩阵单元格将包含1，否则为0。

示例

   A        T
00000001 01111101 
01010100 01100101 
10010111 00010100 
10110000 00011000 <-- This matrix is transposed
11000100 00111110 
10000011 10101111 
11110101 11000100 
10100000 01100010 

二进制和传统乘法的结果：

 Binary  Traditional
11000100  11000100
11111111  32212121
11111111  32213421
11111111  21112211
11101111  22101231
11001111  11001311
11111111  54213432
11001111  11001211

问题

如何以最有效的方式按上述描述相乘这些矩阵？

P.S

我尝试利用二进制 and（即 & 运算符）而不是对单独位进行乘法运算，在这种情况下，我必须为乘法准备数据：

ulong u;

u = T & 0xFF;
u = (u << 00) + (u << 08) + (u << 16) + (u << 24)
  + (u << 32) + (u << 40) + (u << 48) + (u << 56);

现在，通过对两个整数A和u执行二进制and操作，将产生以下结果：

   A        u        R        C
00000001 01111101 00000001    1
01010100 01111101 01010100    3
10010111 01111101 00010101    3
10110000 01111101 00110000    2
11000100 01111101 01000100    2
10000011 01111101 00000001    1
11110101 01111101 01110101    5
10100000 01111101 00100000    1

在上面的示例中，R 包含 A 位乘以 u 位的结果，为了得到最终值，我们必须对一行中的所有位进行求和。请注意，列 C 包含与上述结果为 Traditional 矩阵乘法的第一列中找到的值相等的值。问题在于，在此步骤中，我必须操作单独的位，我认为这是次优的方法，我已经通过http://graphics.stanford.edu/~seander/bithacks.html 查找了一种并行执行的方法，但没有成功，如果有任何人有关于如何将位于 R 列中的值“展平”和“合并”到结果 64 位矩阵中的任何想法，我将不胜感激，期待你的回复。

感谢您，

编辑

非常感谢David Eisenstat，最终算法如下：

var A = ...;
var T = ...; // T == transpose(t), t is original matrix, algorithm works with transposed matrix

var D = 0x8040201008040201UL;

U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D); T = (T << 8) | (T >> 56); D = (D << 8) | (D >> 56);
U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & D);

以下代码段：

    public static void Main (string[] args){
        ulong U;
        var Random = new Xor128 ();

        var timer = DateTime.Now;

        var A = Random.As<IUniformRandom<UInt64>>().Evaluate();
        var T = Random.As<IUniformRandom<UInt64>>().Evaluate();

        var steps = 10000000;

        for (var i = 0; i < steps; i++) {
            ulong r = 0;

            var d = 0x8040201008040201UL;

            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d); T = (T << 8) | (T >> 56); d = (d << 8) | (d >> 56);
            U = A & T; U |= U >> 1; U |= U >> 2; U |= U >> 4; U &= 0x0101010101010101UL; U = (U << 8) - U; r |= (U & d);
        }

        Console.WriteLine (DateTime.Now - timer);


        var m1 = new Int32[8,8];
        var m2 = new Int32[8,8];
        var m3 = new Int32[8,8];

        for (int row = 0; row < 8; row++) {
            for (int col = 0; col < 8; col++) {
                m1 [row, col] = Random.As<IUniformRandom<Int32>> ().Evaluate(0, 1);
                m2 [row, col] = Random.As<IUniformRandom<Int32>> ().Evaluate(0, 1);
                m3 [row, col] = Random.As<IUniformRandom<Int32>> ().Evaluate(0, 1);
            }
        }

        timer = DateTime.Now;

        for (int i = 0; i < steps; i++) {
            for (int row = 0; row < 8; row++) {
                for (int col = 0; col < 8; col++) {
                    var sum = 0;

                    for (int temp = 0; temp < 8; temp++) {
                        sum += m1 [row, temp] * m2 [temp, row];
                    }

                    m3 [row, col] = sum;
                }
            }
        }

        Console.WriteLine (DateTime.Now - timer);
    }

展示以下结果给我：

00:00:02.4035870
00:00:57.5147150

感谢大家，这意味着在Mac OS X/Mono下的性能提升了23倍。

- Lu4

---

也许你应该只标记一种语言... - PlasmaHH

嗯，我不确定我理解你的意思，对我来说两个矩阵只是两个独立的随机数，这样可以吗？ - Lu4

如果A是你的矩阵，T是A的转置，则A和T绝对不是独立的。 - evenex_code

不是A的转置，而是T的转置。T是另一个矩阵，它被转置以使其适用于乘法，在矩阵乘法中，将第一个矩阵的行与第二个矩阵的列进行交叉相乘，如果它被转置，则可以将第一个矩阵的行与第二个矩阵的行进行交叉相乘，这种技巧允许使用二进制位操作进行矩阵乘法... - Lu4

我觉得有必要澄清这一点，因为我和ryyker都不确定你的意思是什么。 - evenex_code

Answer 1

我不确定这是最有效的方法，但以下是一个可以尝试的方法。下面这个指令序列计算了矩阵 A * T' 的主对角线。将矩阵 T 和向量 D 都循环右移 8 位，然后重复执行该过程 7 次。

// uint64_t A, T;
uint64_t D = UINT64_C(0x8040201008040201);
uint64_t P = A & T;
// test whether each byte is nonzero
P |= P >> 1;
P |= P >> 2;
P |= P >> 4;
P &= UINT64_C(0x0101010101010101);
// fill each nonzero byte with ones
P *= 255;  // or P = (P << 8) - P;
// leave only the current diagonal
P &= D;

- David Eisenstat

---

看起来非常优化，我之前的做法远不如此优化。我认为http://graphics.stanford.edu/~seander/bithacks.html缺少你的解决方案，非常感谢！ - Lu4

Answer 2

如果您正在寻找一种并行进行密集矩阵乘法的方法，请将结果矩阵分成块，并同时计算每个块。http://en.wikipedia.org/wiki/Block_matrix#Block_matrix_multiplication

- evenex_code

Answer 3

不清楚您正在使用什么数据结构、哪种语言（是的，我知道你说“任何语言”），以及您要优化什么（速度？内存？等等）。所有这些都可能对您的解决方案产生深远的影响。

以下是一些示例： 1. 如果使用的是C/C++，并且矩阵在内存中是连续的位。每行/列映射到UINT8。在这种情况下，将一行与一列相乘可以简化为执行8位按位 &，并检查结果是否大于0（不需要求和位）。这只需要2个处理器指令。 2. 如果必须进行逐位操作，请使用按位“或”（|）而不是+。有些语言可能会进行惰性计算，在遇到第一个“1”时停止。 3. 如果可以多线程，则可以加快计算速度。

顺便说一下，我假设您需要处理大量矩阵，否则我将使用直接、易读的代码。我的猜测是，即使有很多矩阵，性能收益也将微不足道。

- bavaza

---

目前我正在使用C#，但是我正在将性能关键部分移植到OpenCL（基本上是C），我不要求在OpenCL中提供算法，但任何C / C ++ / C＃或Java的解决方案都会很有用。我正在使用64位无符号整数进行工作。对于我来说，具有性能最优算法至关重要。即使从矩阵乘法中删掉一个多余的指令也可能导致巨大的性能提升，因为目前我每秒执行约500亿次矩阵乘法，目前我正在按位操作，如“P.S.”中所述，这就是基本情况，谢谢。 - Lu4

Answer 4

如果你允许比C/C++更低级别的构建，则SSE/AVX机器指令与内置编译器函数一起可以编写更快的代码（根据我进行的某些基准测试，速度提高了4倍）。你需要使用非标准向量变量（至少由GCC、ICC、CLang支持）：

using epu = uint8_t __attribute__((vector_size(16)));

我正在使用类似于以下的类：

class BMat8 {
    [...]
  private:
    uint64_t _data;
};

然后，下面的代码应该可以实现你想要的功能。

static constexpr epu rothigh { 0, 1, 2, 3, 4, 5, 6, 7,15, 8, 9,10,11,12,13,14};
static constexpr epu rot2    { 6, 7, 0, 1, 2, 3, 4, 5,14,15, 8, 9,10,11,12,13};

inline BMat8 operator*(BMat8 const& tr) const {
  epu x = _mm_set_epi64x(_data, _data);
  epu y = _mm_shuffle_epi8(_mm_set_epi64x(tr._data, tr._data), rothigh);
  epu data {};
  epu diag =  {0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80,
               0x80,0x01,0x02,0x04,0x08,0x10,0x20,0x40};
  for (int i = 0; i < 4; ++i) {
    data |= ((x & y) != epu {}) & diag;
    y    = _mm_shuffle_epi8(y, rot2);
    diag = _mm_shuffle_epi8(diag, rot2);
  }
  return BMat8(_mm_extract_epi64(data, 0) | _mm_extract_epi64(data, 1));
}

特别地，使用128位寄存器，我能够同时执行两次迭代。

- hivert

Answer 5

在x86-64上，可以通过我在此处描述的解决方案相当高效地实现严格的布尔代数。

https://dev59.com/ZW3Xa4cB1Zd3GeqPisz8#55307540

唯一的区别在于，转置矩阵的数据也需要按列提取并重新组合成行，然后再进行每个64位乘积。幸运的是，可以使用BMI2指令进行并行位提取，并通过GCC中的内在函数_pext_u64轻松实现此操作。

uint64_t mul8x8T (uint64_t A, uint64_t B) {

    const uint64_t COL = 0x0101010101010101;

    uint64_t C = 0;

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


uint64_t torow (uint64_t c) {
    const uint64_t ROW = 0x00000000000000FF; // mask of the first row
    const uint64_t COL = 0x0101010101010101; // mask of the first column

    // select bits of c in positions marked by COL,
    // and pack them consecutively
    // last 'and' is included for clarity and is not 
    // really necessary 
    return _pext_u64(c, COL) & ROW;
}

在不支持这个特定指令的处理器上，一种可能的解决方案是调整通常用于打包的典型位操作技巧，例如在使用64位乘法进行字节位逆序时使用的经典位顺序翻转。

https://graphics.stanford.edu/~seander/bithacks.html#ReverseByteWith64BitsDiv

使用掩码和整数乘法与某些常量相乘，得到一个四字节的结果，其中包含打包结果作为位子串，可以使用位移和掩码提取。

思路是将乘法步骤视为并行位移，其中输入中的每个位都被不同的量移动，该量由常量指定。只要两个数字的跨度在结果的某个位置上不发生碰撞，即只要乘法的每个部分和更新结果的不同位位置，这就总是可能的。这避免了任何潜在的进位，使逐位求和等价于位并行 OR（或 XOR）。

uint64_t torow (uint64_t c) {
    const uint64_t ROW = 0x00000000000000FF; // select 8 lowest consecutive bits to get the first row
    const uint64_t COL = 0x0101010101010101; // select every 8th bit to get the first column
    const uint64_t DIA = 0x8040201008040201; // select every 8+1 bit to obtain a diagonal

    c *= ROW; // "copies" first column to the rest
    c &= DIA; // only use diagonal bits or else there will be position collisions and unexpected carries
    c *= COL; // "scatters" every bit to all rows after itself; the last row will now contain the packed bits
    return c >> 56; // move last row to first & discard the rest
}

这个函数还可以使用更低强度的操作进行其他可能的替代实现，最快的实现取决于目标架构。

- paperclip optimizer