algorithm - 二进制矩阵乘法位旋转破解

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

＃问题

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

#PS

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

ulong u;

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

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

   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相乘的结果，为了获得最终值，我们必须将所有位排成一行。请注意，该列包含的值等于上面结果矩阵乘法的第一列中的值。问题是，在这一步中，我必须对我认为次优方法的单独位进行操作，我已通读http://graphics.stanford.edu/~seander/bithacks.html寻找一种方法并行执行此操作但没有运气，如果有人对如何将列中的值“展平”和“合并”成生成的 64 位矩阵有任何想法，如果您给我写几行，我将不胜感激，AusumCTraditionalR

谢谢，

＃编辑

非常感谢 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 倍性能提升，谢谢大家

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;

Answer 2

目前尚不清楚您使用的是哪种数据结构、哪种语言（是的，我知道您说“任何语言”）以及您要优化的内容（速度？内存？）等。所有这些都可能对您的解决方案。

一些例子：

• 假设这是 C/C++，你的矩阵是内存中的连续位。每行/列映射到一个 UINT8。在这种情况下，将行与列相乘减少为执行 8 位按位 &，并检查结果是否大于 0（无需对位求和）。这需要 2 个处理器指令。
• 如果您被迫进行逐位操作，请使用按位“或”( |) 而不是+. 有些语言可能会懒惰地评估这一点，在遇到的第一个“1”处停止。
• 如果您可以多线程，则可以加快计算速度。

顺便说一句，我假设您有很多矩阵要处理，否则我会使用直接且可读的代码。我的猜测是，即使有很多矩阵，性能的提升也可以忽略不计。

Answer 3

如果您正在寻找一种并行执行密集矩阵乘法的方法，请将结果矩阵划分为块并并行计算每个块。

http://en.wikipedia.org/wiki/Block_matrix#Block_matrix_multiplication

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 位寄存器，我可以一次进行两次迭代。

Answer 5

使用我在此处描述的解决方案，可以在 x86-64 上非常有效地实现严格布尔代数的解决方案：

https://stackoverflow.com/a/55307540/11147804

唯一的区别是转置矩阵中的数据也需要按列提取并在每个 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
}

此功能还有其他可能的替代实现，使用更多强度较低的操作，其中最快的将取决于目标架构。