你将如何转置二进制矩阵

chenzoff

已于 2022-08-16 10:27:10 修改

阅读量396

点赞数 1

文章标签：矩阵线性代数

于 2022-08-16 10:16:57 首次发布

原文链接：https://www.likecs.com/ask-1086029.html

版权

【问题标题】：你将如何转置二进制矩阵？(How would you transpose a binary matrix?)
【发布时间】：2015-10-22 21:57:51
【问题描述】：

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

例如下面的矩阵：

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,
};

注意：我更喜欢一种算法，它可以使用任意大小的矩阵进行计算，但我也对只能处理特定大小的算法感兴趣。

【问题讨论】：

我不明白转置后的输出：为什么第一行是0b00000000 而不是0b00000001？为什么第二行是0b00000100 而不是0b00000010？ ...
我看不出如何才能真正避免手动计算每一位。结果的每一行都来自源代码的每一行。这确实会阻止任何有用的并行性......
@m.s.因为第一列全是零（矩阵是 3x7 但用 3x8 表示）
由于是3×8矩阵，转置的输出是8×3矩阵。转置意味着列变成行。
从给定的示例中，我看到您总是转置 8x8 矩阵（对于较小的矩阵，您只需用零填充剩余的行）。在这种情况下，如果您的代码将在 64 位 CPU 上运行，则存在一个非常有效的算法。它在 Knuth 的“计算机编程的艺术”卷中有所描述。 4a 第 7.1.3 章。您可以在此页面上的函数“flipDiagA1H8”中找到一个实现：chessprogramming.wikispaces.com/Flipping+Mirroring+and+Rotating

标签： c++ math matrix binary transpose

【解决方案1】：

我花了更多时间寻找解决方案，并且找到了一些好的解决方案。

SSE2 方式

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

这个解决方案的灵感来自this blog post by mischasan，并且远远优于我迄今为止对这个问题提出的所有建议。

这个想法很简单：

#include <emmintrin.h>
将 16 个 uint8_t 变量打包成一个 __m128i
使用_mm_movemask_epi8 获取每个字节的MSB，生成uint16_t
使用_mm_slli_epi64将128位寄存器移一位
重复直到你得到所有 8 个uint16_ts

通用 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_ts 打包到一个 32 位变量中
屏蔽每个字节的第一位（使用0x80808080）
乘以0x02040810
取乘法的高 32 位的 4 个 LSB
通常，您可以屏蔽每个字节中的第 N 位（将掩码右移 N 位）并乘以幻数，左移 N 位。这里的好处是，如果您的编译器足够聪明，可以展开循环，那么掩码和“幻数”都会成为编译时常量，因此移动它们不会导致任何性能损失。最后一个 4 位序列有一些问题，因为丢失了一个 LSB，所以在这种情况下，我需要将输入左移 8 位，并使用与第一个 4 位序列相同的方法。

如果您使用两个 4×8 块进行此操作，那么您可以完成一个 8x8 块并排列生成的位，以便所有内容都放在正确的位置。

【问题讨论】：

@étale-cohomology 是的，见github.com/Venemo/fecmagic/blob/master/src/binarymatrix.h
@étale-cohomology 从长远来看，SSE 版本不值得付出努力，因为编译器可以优化“通用”版本以比我使用 SSE 编写的版本更快。 :)

【解决方案2】：

我的建议是，不要进行转置，而是在矩阵数据中添加一位信息，指示矩阵是否转置。

现在，如果要将转置矩阵与向量相乘，则与将左侧的矩阵乘以向量（然后转置）相同。这很简单：只需对 8 位数字进行一些 xor 操作即可。

然而，这使得其他一些操作变得复杂（例如添加两个矩阵）。但是在评论中你说乘法正是你想要优化的。

【问题讨论】：

抱歉，这不是一个可接受的解决方案。我需要实际转置矩阵，因为我需要在交织器的输出端使用它。 :) 无论如何，我欢迎一个关于如何在不实际转置矩阵的情况下执行向量乘法的示例。
这很聪明（因此 +1），但是我的用例需要计算 A*T(A) ;/

【解决方案3】：

这是 Jay Foad 给我的关于快速布尔矩阵的电子邮件文本转置：

布尔转置算法的核心是一个我称之为transpose8x8 的函数，它转置封装在一个 64 位字中的 8x8 布尔矩阵（按照从 MSB 到 LSB 的行主要顺序）。要转置宽度和高度为 8 的倍数的任何矩形矩阵，请将其分解为 8x8 块，单独转置每个块并将它们存储在输出中的适当位置。要加载 8x8 块，您必须加载 8 个单独的字节并将它们移位和 OR 成一个 64 位字。存储也是如此。

transpose8x8 的普通 C 实现依赖于这样一个事实：任何平行于前导对角线的对角线上的所有位都向上/向下和向左/向右移动相同的距离。例如，前导对角线上方的所有位都必须向左移动一位，向下移动一位，即在打包的 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 指令，您可以实现完美的 shuffle，并使用它来执行 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 个单独的字节。）

【问题讨论】：

什么是PDEP 和PEXT 指令？什么是力量？您能否发布您的完整实现，以便我可以针对我的实现进行基准测试？
PDEP 和 PEXT 是
@Venemo POWER 是现代 IBM 大型机（除了 Z/Architecture）上使用的 ISA。它基于 PowerPC。此指令不适用于 x86 或 ARM 计算机。

【解决方案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 次！我们可以单独对每个字节进行操作。

问题

正如 cmets 中有人指出的那样：Translate(Translate(M)) != M 所以如果你需要这个，你需要做一些额外的工作。

可以通过直接映射 uint64_t 而不是 uint8_t[8] 数组来提高性能，因为它省略了“安全复制”以防止对齐问题。

【问题讨论】：

查找表有多大？我的意思是，在你的代码中它只有 256 长，但不确定这个数字来自哪里。
查找表为 2kiB（2048 字节）。它有 256 个 64 位（8 字节）的条目。原因是您使用的是 uint8_t，它是 8 位 --> 256 个可能的值 (2^8)，一个 uint8_t 可以包含 0 到 255。此代码的限制是它最多只支持 8x8大小。
这个限制也是它能够工作的原因，一个 uint64_t 有 64 位，方便地等于 8x8，这使我们能够使用一个 uint64_t 作为一个 8x8 位矩阵并执行一些数学运算和或对其进行位操作。

【解决方案5】：

我添加了一个新的遮阳篷，而不是编辑我原来的遮阳篷，以使其更加明显（遗憾的是没有评论权）。

在您自己的 awnser 中添加第一个中没有的附加要求：它必须在 ARM Cortex-M 上工作

我确实在我原来的 awnser 中为 ARM 提出了一个替代解决方案，但由于它不是问题的一部分并且似乎离题（主要是因为 C++ 标签）而忽略了它。

ARM 专用解决方案 Cortex-M：

部分或大部分 Cortex-M 3/4 有一个位带区域，可以完全满足您的需要，它将位扩展为 32 位字段，该区域可用于执行原子位操作。

如果您将阵列放在位带区域中，它将在位带区域中有一个“爆炸”镜像，您可以在其中对位本身使用移动操作。如果您创建一个循环，编译器肯定能够展开和优化以仅移动操作。

如果你真的想要，你甚至可以设置一个 DMA 控制器来处理一整批转置操作，并完全从 cpu 中卸载 :)

也许这仍然对你有帮助。

【问题讨论】：

感谢您的帖子！我不知道位带。老实说，我不确定它是否会加快我的方法。一次处理四个压缩的 8 位数字似乎优于处理单个位。
无论如何，你的两个答案都得到了我的支持，我很感激你花时间思考我的问题。 :)
是的，一次在多个位上操作的转置当然会更快。我确实想提一下这个选项，因为它可能对矩阵乘法、位提取/测试/设置等其他部分有用，因为它可以直接对位进行操作，因此您可以保存所有移位和屏蔽操作，它只会打开很多其他可能性作为我在示例中使用 DMA 进行大批量的示例:)

【解决方案6】：

这有点晚了，但我今天偶然发现了这个交换。如果您查看 Hacker's Delight，第 2 版，从第 141 页开始，有几种算法可以有效地转置布尔数组。

它们非常有效：我的一位同事获得了大约 10 倍的系数在 X86 上与幼稚编码相比，速度更快。

【问题讨论】：

您能说得更具体些吗？通常在 StackOverflow 上，您要么给出完整答案，要么给出完整答案的链接。
书里的答案占了好几页，所以我舍不得
...不愿意尝试正确输入，更不用说可能侵犯作者的版权。
至少尝试描述一下该算法的基本思想，如何以及为什么它比当前接受的答案更好，等等。
Jay Foad (www.dyalog.com) 好心让我发一封电子邮件

【解决方案7】：

这是我在 gitub 上发布的内容 (mischasan/sse2/ssebmx.src) 将 INP() 和 OUT() 更改为使用归纳变量可以分别保存一个 IMUL。 AVX256 的速度是原来的两倍。 AVX512 不是一个选项，因为没有 _mm512_movemask_epi8()。

#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);
}

HTH。

【问题讨论】：

【解决方案8】：

受 Roberts 回答的启发，Arm Neon 中的多项式乘法可用于分散比特 --

inline poly8x16_t mull_lo(poly8x16_t a) {
     auto b = vget_low_p8(a);
     return vreinterpretq_p8_p16(vmull_p8(b,b));
}
inline poly8x16_t mull_hi(poly8x16_t a) {
     auto b = vget_high_p8(a);
     return vreinterpretq_p8_p16(vmull_p8(b,b));
}
auto a = mull_lo(word);
auto b = mull_lo(a), c = mull_hi(a);
auto d = mull_lo(b), e = mull_hi(b);
auto f = mull_lo(c), g = mull_hi(c);

那么vsli可以用来成对组合位。

auto ab = vsli_p8(vget_high_p8(d), vget_low_p8(d), 1);
auto cd = vsli_p8(vget_high_p8(e), vget_low_p8(e), 1);
auto ef = vsli_p8(vget_high_p8(f), vget_low_p8(f), 1);
auto gh = vsli_p8(vget_high_p8(g), vget_low_p8(g), 1);
auto abcd = vsli_p8(ab, cd, 2);
auto efgh = vsli_p8(ef, gh, 2);
return vsli_p8(abcd, efgh, 4);

Clang 优化此代码以避免 vmull2 指令，大量使用 ext q0,q0,8 到 vget_high_p8。

迭代方法可能不仅速度更快，而且使用的寄存器更少，并且可以简化 2 倍或更多的吞吐量。

// 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

【问题讨论】：