如何转置二进制矩阵?

11

我在C++中有二进制矩阵,使用8位值的向量表示。

例如,以下矩阵:

I have binary matrices in C++ that I repesent with a vector of 8-bit values.

For example, the following matrix:

1 0 1 0 1 0 1
0 1 1 0 0 1 1
0 0 0 1 1 1 1

表示为:

const uint8_t matrix[] = {
    0b01010101,
    0b00110011,
    0b00001111,
};

我这样做的原因是,这样计算这种矩阵和一个8位向量的乘积变得非常简单和高效(每行只需进行一次按位与运算和奇偶校验计算),这比逐个计算每个比特位要好得多。

我现在正在寻找一种高效的方法来转置这样的矩阵,但是我还没有找到如何在不手动计算每个比特位的情况下完成它。

仅供澄清,对于上述示例,我希望从转置中获得以下结果:

const uint8_t transposed[] = {
    0b00000000,
    0b00000100,
    0b00000010,
    0b00000110,
    0b00000001,
    0b00000101,
    0b00000011,
    0b00000111,
};

注意:我更希望有一种可以计算任意大小矩阵的算法,但我也对那些只能处理特定大小矩阵的算法感兴趣。


1
我不理解转置输出:为什么第一行是0b00000000而不是0b00000001?为什么第二行是0b00000100而不是0b00000010?... - m.s.
2
我看不出你如何真正避免手动计算每一位。你结果的每一行都有来自源数据的每一行的一位。这确实防止了任何有用的并行处理... - The Archetypal Paul
由于这是一个3×8矩阵,所以转置的输出是一个8×3矩阵。转置意味着列变成了行。 - Venemo
你需要实际转置数据吗?这样非常低效。你能否将数据包装在接口中,并只标记数据已转置?这可以帮助优化一些矩阵操作。 - Conor
4
从给定的示例中我看到,您总是转置8x8矩阵(对于较小的矩阵,您只需用零填充剩余的行)。在这种情况下,如果您的代码能够在64位CPU上工作,则存在一种相当高效的算法。它在Knuth的《计算机程序设计艺术》第4a卷7.1.3章节中有描述。您可以在此页面的“flipDiagA1H8”函数中找到一个实现:https://chessprogramming.wikispaces.com/Flipping+Mirroring+and+Rotating - Evgeny Kluev
显示剩余9条评论
8个回答

8

我花了更多的时间寻找解决方案,终于找到了一些好的方法。

SSE2方法

在现代x86 CPU上,可以使用SSE2指令非常高效地转置二进制矩阵。使用这样的指令,可以处理一个16×8的矩阵。

这个解决方案受到mischasan的博客文章的启发,比我到目前为止得到的每个建议都要好得多。

思路很简单:

  • #include <emmintrin.h>
  • 将16个uint8_t变量打包成一个__m128i
  • 使用_mm_movemask_epi8获取每个字节的MSB,生成一个uint16_t
  • 使用_mm_slli_epi64将128位寄存器向左移动一位
  • 重复以上步骤直到获得所有8个uint16_t

通用32位解决方案

不幸的是,我还需要让它在ARM上工作。在实现SSE2版本后,只需找到NEON的等效版本即可,但是Cortex-M CPU(与Cortex-A相反)没有SIMD功能,因此目前对我来说NEON并不太有用。

注意:由于Cortex-M没有本地64位算术,因此我无法使用任何建议中将8x8块视为uint64_t的方法。大多数具有Cortex-M CPU的微控制器也没有太多内存,因此我更喜欢在不使用查找表的情况下完成所有这些操作。

经过一些思考,可以使用普通的32位算术和一些巧妙的编码来实现相同的算法。这样,我可以一次处理4×8个块。这是由一位同事提出的建议,其中魔法在于32位乘法的工作方式:您可以找到一个32位数字,然后每个字节的MSB都会出现在结果的上32位中。

  • 将4个uint8_t打包在一个32位变量中
  • 屏蔽每个字节的第1位(使用0x80808080进行屏蔽)
  • 乘以0x02040810
  • 取乘法的上32位的4个LSB(最低有效位)
  • 通常,您可以屏蔽每个字节中的第N位(将掩码右移N位),然后乘以魔术数,并将其左移N位。这里的优点是,如果您的编译器足够聪明,可以展开循环,则掩码和“魔术数”都成为编译时常量,因此移动它们不会产生任何性能损失。对于最后一组4位,存在一些问题,因为这样会丢失一个LSB,因此在这种情况下,我需要将输入向左移动8位,并使用与第一组4位相同的方法。

如果使用两个4×8块,则可以得到一个8x8块,并安排结果位,使所有内容都落在正确的位置。


1
@étale-cohomology 是的,请参见 https://github.com/Venemo/fecmagic/blob/master/src/binarymatrix.h - Venemo
1
最终结果表明,SSE版本并不值得花费精力,因为编译器可以将“通用”版本优化为比我使用SSE写的更快。 :) - Venemo

6
以下是Jay Foad给我的邮件内容,有关快速布尔矩阵转置:
布尔转置算法的核心是一个我称之为transpose8x8的函数,它可以转置一个8x8的布尔矩阵,该矩阵以64位字(从高位到低位按行优先顺序)打包。要转置任何宽度和高度均为8的倍数的矩阵,请将其分解为8x8块,逐个转置并将它们存储在输出中的适当位置。要加载8x8块,您必须加载8个单独的字节,并将它们移位和OR运算到一个64位字中。存储也需要相同的操作。
使用普通的C实现transpose8x8 依赖于这样一个事实:平行于主对角线的任何对角线上的所有比特都向上/向下和左/右移动相同的距离。例如,紧靠着主对角线上方的所有比特都要向左移动一位并向下移动一位,即在打包的64位字中向右移动7位。这导致了以下算法:
transpose8x8(word) {

  return
    (word & 0x0100000000000000) >> 49 // top right corner

  | (word & 0x0201000000000000) >> 42

  | ...

  | (word & 0x4020100804020100) >> 7 // just above diagonal

  | (word & 0x8040201008040201) // leading diagonal

  | (word & 0x0080402010080402) << 7 // just below diagonal

  | ...
  | (word & 0x0000000000008040) << 42

  | (word & 0x0000000000000080) << 49; // bottom left corner

}

这个实现比以前的实现快了大约10倍,以前的实现是从内存中逐个复制源字节中的每个位并将其合并到内存中的目标字节中。

或者,如果您有PDEP和PEXT指令,可以实现完美洗牌,并使用它来执行转置,如Hacker's Delight中所述。 这样会快得多(但我手头没有时间):

shuffle(word) {
    return pdep(word >> 32, 0xaaaaaaaaaaaaaaaa) | pdep(word, 0x5555555555555555);
} // outer perfect shuffle

transpose8x8(word) { return shuffle(shuffle(shuffle(word))); }

POWER的vgbbd指令可以在单个指令中有效地实现整个transpose8x8函数(由于它是一个128位矢量指令,因此它会分别在低64位和高64位上执行两次)。与普通的C实现相比,这使速度提高了约15%。(只有15%是因为虽然位操作更快,但总运行时间现在主要受到加载8字节并将其组装成transpose8x8参数的时间以及获取结果并将其存储为8个单独字节的时间的影响。)


1
PDEPPEXT指令是什么?POWER是什么?您能否发布您的完整实现,以便我可以对我的实现进行基准测试? - Venemo
PDEP和PEXT是 - Robert Bernecky
@Venemo POWER是现代IBM大型机上使用的ISA(除了Z/Architecture)。它基于PowerPC。这个指令在x86或ARM计算机上不可用。 - fuz

5
我的建议是,不要对矩阵进行转置,而是向您的矩阵数据添加一位信息,指示矩阵是否被转置。
现在,如果您想将一个转置矩阵与一个向量相乘,它将等同于在向量左边乘以矩阵(然后转置)。这很容易:只需要对您的8位数字进行一些xor操作。
然而,这使得其他一些操作变得复杂(例如,将两个矩阵相加)。但是,在您的评论中,您说乘法正是您想要优化的。

1
抱歉,但这不是一个可接受的解决方案。我需要实际转置矩阵,因为我需要它在交织器的输出处。 :) 无论如何,我很欢迎提供一个示例,演示如何执行向量乘法而不实际转置矩阵。 - Venemo
这很聪明(因此+1),但是我的用例需要计算A*T(A) ;/ - Brian Vandenberg

4
我的建议是使用查找表来加速处理。
另一件要注意的事情是,根据当前矩阵定义,最大尺寸将为8x8位。这适用于uint64_t,因此我们可以利用这一点,特别是在使用64位平台时。
我已经设计出一个简单的示例,使用查找表,您可以在下面找到并使用:http://www.tutorialspoint.com/compile_cpp11_online.php在线编译器。
示例代码
#include <iostream>
#include <bitset>
#include <stdint.h>
#include <assert.h>

using std::cout;
using std::endl;
using std::bitset;

/* Static lookup table */
static uint64_t lut[256];

/* Helper function to print array */
template<int N>
void print_arr(const uint8_t (&arr)[N]){
    for(int i=0; i < N; ++i){
        cout << bitset<8>(arr[i]) << endl;
    }
}

/* Transpose function */

template<int N>
void transpose_bitmatrix(const uint8_t (&matrix)[N], uint8_t (&transposed)[8]){
    assert(N <= 8);

    uint64_t value = 0;
    for(int i=0; i < N; ++i){
        value = (value << 1) + lut[matrix[i]];
    }

    /* Ensure safe copy to prevent misalignment issues */
    /* Can be removed if input array can be treated as uint64_t directly */
    for(int i=0; i < 8; ++i){
        transposed[i] = (value >> (i * 8)) & 0xFF;
    }
}

/* Calculate lookup table */
void calculate_lut(void){
    /* For all byte values */
    for(uint64_t i = 0; i < 256; ++i){
        auto b = std::bitset<8>(i);
        auto v = std::bitset<64>(0);

        /* For all bits in current byte */
        for(int bit=0; bit < 8; ++bit){
            if(b.test(bit)){
                v.set((7 - bit) * 8);
            }
        }

        lut[i] = v.to_ullong();
    }
}

int main()
{
    calculate_lut();

    const uint8_t matrix[] = {
        0b01010101,
        0b00110011,
        0b00001111,
    };

    uint8_t transposed[8];

    transpose_bitmatrix(matrix, transposed);
    print_arr(transposed);

   return 0;
}

工作原理

您的3x8矩阵将被转置为8x3矩阵,并表示为一个8x8数组。问题在于,您想要将位(“水平”表示)转换为垂直表示,并分散在几个字节中。

正如我上面所提到的,我们可以利用输出(8x8)总是适合于uint64_t的事实。我们将利用此功能,因为现在我们可以使用uint64_t来写入8字节数组,但我们也可以用它来添加、异或等,因为我们可以对64位整数执行基本算术运算。

您3x8矩阵(输入)中的每个条目都有8位宽度,为了优化处理,我们首先生成256条目查找表(每个字节值)。该条目本身是一个uint64_t,将包含位的旋转版本。

例如:

字节= 0b01001111 = 0x4F
lut [0x4F] = 0x0001000001010101 =(uint8_t []){0,1,0,0,1,1,1,1}

现在进行计算:

对于计算,我们使用uint64_t,但请记住,在水下它将表示uint8_t [8]数组。我们简单地移动当前值(从0开始),查找第一个字节并将其添加到当前值中。

这里的“神奇”之处在于,查找表中uint64_t的每个字节只会是1或0,因此它仅设置最低有效位(每个字节)。移动uint64_t将移动每个字节,只要我们确保不要这样做超过8次!,我们可以单独对每个字节执行操作。

问题

正如某人在评论中指出的那样:Translate(Translate(M))!= M,所以如果您需要这个,则需要进行一些额外的工作。

通过直接映射uint64_t而不是uint8_t [8]数组,性能可以得到改进,因为它省略了一个“安全复制”,以防止对齐问题。


查找表有多大?我的意思是,在你的代码中它只有256个长度,但不确定这个数字从哪里来。 - Venemo
查找表的大小为2kiB(2048字节)。它有256个64位(8字节)的条目。原因是您正在使用uint8_t,它是8位-> 256个可能值(2 ^ 8),一个uint8_t可以包含0到255。此代码的限制是它仅支持最大尺寸为8x8。 - AlfaVector
限制也是这能够工作的原因,uint64_t有64位,恰巧等于8x8,使我们可以将uint64_t用作8x8位矩阵,并对其执行一些数学或位运算。 - AlfaVector

3

我添加了一个新的答案,而不是编辑我的原始答案,以使其更加可见(不幸的是没有评论权限)。

在您自己的答案中,您添加了第一个答案中不存在的附加要求:它必须在ARM Cortex-M上运行

在我的原始答案中,我提出了一种针对ARM的替代解决方案,但由于它不属于问题范围,并且似乎与C ++标签无关,因此省略了它。

ARM特定解决方案Cortex-M:

一些或大多数Cortex-M 3/4具有位带区域,可以用于完全满足您的需求,它将位扩展为32位字段,该区域可用于执行原子位操作。

如果您将数组放入位带区域中,则会在位带区域中生成一个“扩展”镜像,您可以直接在位本身上进行移动操作。如果您制作循环,编译器肯定能够展开和优化到仅使用移动操作。

如果您真的想要,甚至可以设置DMA控制器来处理一整批转置操作,并且通过一点努力完全卸载它,从而完全脱离CPU :)

也许这仍然可以帮助您。


谢谢您的帖子!我之前不知道位带操作。老实说,我不确定它是否能加速我的方法。一次处理四个打包的8位数字似乎比单个位更优越。 - Venemo
您IP地址为143.198.54.68,由于运营成本限制,当前对于免费用户的使用频率限制为每个IP每72小时10次对话,如需解除限制,请点击左下角设置图标按钮(手机用户先点击左上角菜单按钮)。 - Venemo
1
当然,对于一次转置多个位的操作会更快。我想提一下这个选项,因为它可能对其他部分有用,比如矩阵乘法、位提取/测试/设置等,因为它可以直接操作位,因此您可以节省所有的移位和掩码操作,并且它只是打开了许多其他可能性,就像我在示例中使用DMA来处理大批量数据一样 :) - AlfaVector

2
这是我在GitHub上发布的内容(mischasan/sse2/ssebmx.src): 把INP()和OUT()改为使用归纳变量可以节省一次IMUL。 AVX256速度是它的两倍。 因为没有_mm512_movemask_epi8(),所以AVX512不是一个选项。
#include <stdint.h>
#include <emmintrin.h>

#define INP(x,y) inp[(x)*ncols/8 + (y)/8]
#define OUT(x,y) out[(y)*nrows/8 + (x)/8]

void ssebmx(char const *inp, char *out, int nrows, int ncols)
{
    int rr, cc, i, h;
    union { __m128i x; uint8_t b[16]; } tmp;

    // Do the main body in [16 x 8] blocks:
    for (rr = 0; rr <= nrows - 16; rr += 16)
        for (cc = 0; cc < ncols; cc += 8) {
            for (i = 0; i < 16; ++i)
                tmp.b[i] = INP(rr + i, cc);
            for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
                *(uint16_t*)&OUT(rr, cc + i) = _mm_movemask_epi8(tmp.x);
        }

    if (rr == nrows) return;

    // The remainder is a row of [8 x 16]* [8 x 8]?

    //  Do the [8 x 16] blocks:
    for (cc = 0; cc <= ncols - 16; cc += 16) {
        for (i = 8; i--;)
            tmp.b[i] = h = *(uint16_t const*)&INP(rr + i, cc),
            tmp.b[i + 8] = h >> 8;
        for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
            OUT(rr, cc + i) = h = _mm_movemask_epi8(tmp.x),
            OUT(rr, cc + i + 8) = h >> 8;
    }

    if (cc == ncols) return;

    //  Do the remaining [8 x 8] block:
    for (i = 8; i--;)
        tmp.b[i] = INP(rr + i, cc);
    for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
        OUT(rr, cc + i) = _mm_movemask_epi8(tmp.x);
}

祝好。


2

虽然有点晚,但我今天刚刚偶然发现了这个交流。 如果你查看《Hacker's Delight,第二版》,就会发现从第141页开始,有几个有效地转置布尔数组的算法。

它们非常高效:我的一个同事在X86上与朴素编码相比,获得了约10倍的加速比。


不太愿意尝试正确地打出所有内容,更不用说潜在的侵犯作者版权了。 - Robert Bernecky
至少尝试描述该算法的基本思想,以及它为什么比当前接受的答案更好等方面。 - Venemo
Jay还说,如果你有PDEP和PEXT,你可以实现完美洗牌,并以更快的速度进行转置。详见Hacker's Delight为什么会这样:shuffle(word) { return pdep(word >> 32, 0xaaaaaaaaaaaaaaaa) | pdep(word, 0x5555555555555555); } // outer perfect shuffle transpose8x8(word) { return shuffle(shuffle(shuffle(word))); } - Robert Bernecky
我建议将此内容添加到您的答案中,而不是在评论部分中。另外,这个方法是否比当前接受的答案表现更好? - Venemo
显示剩余4条评论

0
受到Robert的回答的启发,可以利用Arm Neon中的多项式乘法来分散位数--
poly8x8_t transpose(poly8x8_t a) {
   for (int i = 0; i < 3; i++) {
     auto q = vmull_p8(a,a);
     auto top = vreinterpret_u8_p16(vget_high_p16(q));
     auto low = vreinterpret_u8_p16(vget_low_p16(q));
     low = vadd_u8(low, low); // shift left by 1
     low = vadd_u8(low, top); // interleave the bits
     a = vreinterpret_p8_u8(low);
   }
   return a;
}

64位处理器的通用方法
// transpose  bits in 2x2 blocks, first 4 rows
//   x = a b|c d|e f|g h      a i|c k|e m|g o   | byte 0
//       i j|k l|m n|o p      b j|d l|f n|h p   | byte 1
//       q r|s t|u v|w x      q A|s C|u E|w G   | byte 2
//       A B|C D|E F|G H      r B|t D|v F|h H   | byte 3 ...
// ----------------------

auto a = (x & 0x00aa00aa00aa00aaull);
auto b = (x & 0x5500550055005500ull);
auto c = (x & 0xaa55aa55aa55aa55ull) | (a << 7) | (b >> 7);

// transpose 2x2 blocks (first 4 rows shown)
//   aa bb cc dd      aa ii cc kk
//   ee ff gg hh   -> ee mm gg oo
//   ii jj kk ll      bb jj dd ll
//   mm nn oo pp      ff nn hh pp

auto d = (c & 0x0000cccc0000ccccull);
auto e = (c & 0x3333000033330000ull);
auto f = (c & 0xcccc3333cccc3333ull) | (d << 14) | (e >> 14);

// Final transpose of 4x4 bit blocks

auto g = (f & 0x00000000f0f0f0f0ull);
auto h = (f & 0x0f0f0f0f00000000ull);
x = (f & 0xf0f0f0f00f0f0f0full) | (g << 28) | (h >> 28);

在ARM中,现在每个步骤都可以由3条指令组成。
auto tmp = vrev16_u8(x);
tmp = vshl_u8(tmp, plus_minus_1); // 0xff01ff01ff01ff01ull
x = vbsl_u8(mask_1, x, tmp);   // 0xaa55aa55aa55aa55ull

tmp = vrev32_u16(x);
tmp = vshl_u16(tmp, plus_minus_2); // 0xfefe0202fefe0202ull
x = vbsl_u8(mask_2, x, tmp);   // 0xcccc3333cccc3333ull

tmp = vrev64_u32(x);
tmp = vshl_u32(tmp, plus_minus_4); // 0xfcfcfcfc04040404ull
x = vbsl_u8(mask_4, x, tmp);   // 0xf0f0f0f00f0f0f0full

在没有NEON和有限的任意立即数的Arm64上,可以使用以下方法:
uint64_t transpose(uint64_t x) {
    auto b = x ^ (x >> 7);
    b &= 0x00aa00aa00aa00aa;
    x ^= b;
    x ^= b << 7;

    b = x ^ (x >> 14);
    b &= 0x0000cccc0000cccc;
    x ^= b;
    x ^= b << 14;

    b = x ^ (x >> 28);
    b &= 0x00000000f0f0f0f0;
    x ^= b;
    x ^= (uint64_t)b << 28;
    return x;
}

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接