使用按位 AND 和 popcount 而不是实际的 int 或 float 乘法的大 (0,1) 矩阵乘法？
Large (0,1) matrix multiplication using bitwise AND and popcount instead of actual int or float multiplies?
对于大型二进制矩阵（10Kx20K）的乘法，我通常做的是将矩阵转换为浮点矩阵并执行浮点矩阵乘法，因为整数矩阵乘法非常慢（）。
不过这一次，我需要执行数十万次这样的乘法运算，平均而言，即使是毫秒级的性能提升对我来说也很重要。
我想要一个 int 或 float 矩阵作为结果，因为乘积可能包含非 0 或 1 的元素。输入矩阵元素全部为 0 或 1，因此它们可以存储为单个位。
在行向量和列向量之间的内积（产生输出矩阵的一个元素）中，乘法简化为按位与。加法仍然是加法，但我们可以使用人口计数函数添加位，而不是单独循环它们。
一些其他 boolean/binary-matrix 函数或位而不是计算它们，产生位矩阵结果，但这不是我需要的。
这是一个示例代码，显示将问题形成为 std::bitset、AND 和 count 运算比矩阵乘法更快。
#include <iostream>
using std::cout; using std::endl;
#include <vector>
using std::vector;
#include <chrono>
#include <Eigen/Dense>
using Eigen::Map; using Eigen::Matrix; using Eigen::MatrixXf;
#include <random>
using std::random_device; using std::mt19937; using std::uniform_int_distribution;
#include <bitset>
using std::bitset;
using std::floor;
const int NROW = 1000;
const int NCOL = 20000;
const float DENSITY = 0.4;
const float DENOMINATOR = 10.0 - (10*DENSITY);
void fill_random(vector<float>& vec) {
random_device rd;
mt19937 eng(rd());
uniform_int_distribution<> distr(0, 10);
int nnz = 0;
for (int i = 0; i < NROW*NCOL; ++i)
vec.push_back(floor(distr(eng)/DENOMINATOR));
}
void matmul(vector<float>& vec){
float *p = vec.data();
MatrixXf A = Eigen::Map<Eigen::Matrix<float, NROW, NCOL, Eigen::RowMajor>>(p);
cout << "Eigen matrix has " << A.rows() << " rows and " << A.cols() << " columns." << endl;
cout << "Total non-zero values : " << A.sum() << endl;
cout << "The density of non-zero values is " << A.sum() * 1.0 / (A.cols()*A.rows()) << endl;
auto start = std::chrono::steady_clock::now();
MatrixXf B = A.transpose() * A;
auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
cout << "Mat mul took " << end << " ms"<< endl;
// Just to make sure the operation is not skipped by compiler
cout << "Eigen coo ";
for (int i=0; i<10; ++i)
cout << B(0,i) << " ";
cout << endl;
}
void bitset_op(vector<float>& vec) {
// yeah it's not a great idea to set size at compile time but have to
vector<bitset<NROW>> col_major(NCOL);
// right, multiple par for isn't a good idea, maybe in a parallel block
// Doing this for simplicity to profile second loop timing
// converting row major float vec to col major bool vec
#pragma omp parallel for
for (int j=0; j < NCOL; ++j) {
for (int i=0; i < NROW; ++i) {
col_major[j].set(i, vec[i*NCOL + j] && 1);
}
}
auto start = std::chrono::steady_clock::now();
vector<int> coo;
coo.assign(NCOL*NCOL, 0);
#pragma omp parallel for
for (int j=0; j < NCOL; ++j) {
for (int k=0; k<NCOL; ++k) {
coo[j*NCOL + k] = (col_major[j]&col_major[k]).count();
}
}
auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
cout << "bitset intersection took " << end << " ms"<< endl;
// Just to make sure the operation is not skipped by compiler
cout << "biset coo ";
for (int i=0; i<10; ++i)
cout << coo[i] << " ";
cout << endl;
}
int main() {
// Saving to float instead of int to speed up matmul
vector<float> vec;
fill_random(vec);
matmul(vec);
bitset_op(vec);
}
运行 这与：
g++ -O3 -fopenmp -march=native -I. -std=c++11 code.cpp -o code
我得到：
Eigen matrix has 1000 rows and 20000 columns.
Total non-zero values : 9.08978e+06
The density of non-zero values is 0.454489
Mat mul took 1849 ms
Eigen coo 458 206 208 201 224 205 204 199 217 210
bitset intersection took 602 ms
biset coo 458 206 208 201 224 205 204 199 217 210
如您所见，matmul 作为一组 bitset 操作比 Eigen 的 float matmul 快大约 3 倍，这是有道理的。
我想强调的是，我需要在 100K（在 HPC 或云中）上执行此操作，平均毫秒级的性能提升会有所不同。
我不受任何特定库、C++ 标准等的约束。所以请随时回答您认为比使用 GPU 的解决方案更快的任何解决方案，因为我不能将它用于数字原因。
我不确定你可以让编译器为你做多少，如果有的话，而无需使用内在函数或 C++ vector-class 包装器手动矢量化（如 Agner Fog's VCL，如果你的项目是许可证与 GPL 兼容）。还有一些 non-GPLed 包装器。
Cache-blocking 矩阵乘法是一门艺术（在这里很重要），如果你可以使用 Eigen 的现有模板，那将是非常好的，但是一个不同的 class，它在整数上使用按位 and，而不是在浮点数上乘法。我不确定这是否可行。
我做了一些搜索，大多数关于二进制矩阵的文献都是关于产生布尔结果的（包括 SO 问题 like this）。向量 inner-product 是用 AND 作为乘法，但用 XOR 或 OR 作为加法，而不是 popcount。也许有一个 search-term 我错过了它描述 "normal" 矩阵恰好是 (0,1) 矩阵，但乘积不是。
由于每一毫秒都很重要，您可能需要手动对其进行矢量化。
并不是说 vector-integer 东西一般都很慢，只是 只是 vector-integer 相乘比较慢，尤其是与 vector-float 相比最近的 x86 硬件上的 FMA（尤其是英特尔，它在 Haswell 及更高版本上的 FP FMA 吞吐量为每时钟 2x 256b 矢量）。
由于您不需要实际与布尔元素相乘，只需一个 AND（每个时钟吞吐量 3 个向量），这对您来说不是问题。 efficiency-gain 从每个向量中获取更多元素应该足以弥补任何额外成本 per-vector.
当然，这假设一个整数 matmul 实现使用所有相同的 cache-blocking 和其他优化作为等效的 FP matmul，如果你不想（或不想），这就是问题所在知道如何）自己编写，找不到可以为您编写的库。
我只是在回答一个问题，即通过最佳实施， 的效率如何。 标题问题的答案是非常明确的是，使用实际乘法是浪费大量时间，尤其是对于 32 位元素。
存储格式选项：
每字节一个元素(0/1):
- 4 倍密度
float（缓存占用空间/内存带宽/每个向量的元素）
- 易于使用字节洗牌进行转置
- 垂直 ADD 很容易，以防万一（例如，在外循环上进行矢量化，并同时处理多行或多列。如果您有数据，可能会很好（避免最后水平求和）交错的方式使这项工作无需额外的改组。）
每字节4个元素，打包到低半字节:
- 4 倍的单独字节密度
- 使用 AVX2
vpshufb 进行 popcount 非常有效。对于 L1D 缓存中的热输入，理论上，您可以 load/AND/accumulate-a-popcount 每个时钟周期（每个内核）吞吐 128 AND-result 个元素。每个时钟 4 fused-domain 微指令饱和 SKL/HSW front-end 发出每个时钟 4 的带宽，并且不会在 3 个矢量 ALU 端口上出现瓶颈，因为其中一个微指令是纯负载。 （另一个负载 micro-fuses 与 vpand）。当 L2 带宽出现瓶颈时（每个周期约 32B 负载），运行s 在每个时钟 64 个元素。见下文。
- 从整数或 packed-bitmap 创建速度较慢（但如果您以交错顺序将位放入向量以实现高效 pack/unpack 到 in-order 字节，而不是强制位到井然有序）。
- 难以转置（可能比满载还差）
压缩位:
- 8 倍单独字节的密度，每个 AVX2 向量 256 个元素。
- 可以从具有
pmovmskb 的向量创建，用于 non-interleaved 存储顺序。 （不过，对于动态创建不是很有用，因为这会将结果放在整数 reg 中，而不是向量中。交错的 bit-order 可能是最好的，尤其是在转置期间解包）。
- 使用 AVX2 进行 popcount 相当有效：mask / shift+mask / 2x
vpshufb。 (9 fused-domain uops (8 vector-ALU uops) 到 AND + 累积 256 个元素的 popcount（来自 2 row/column 向量），对比 8 uops (6 vector-ALU uops) 4-per-byte 策略（来自 4 row/column 向量）。ALU 端口瓶颈将其限制为来自 L1D 或 L2 的每个时钟 96 个元素。因此，当它在 L2 带宽上出现瓶颈时，这大约是 pack4 策略的 inner-product 吞吐量的 1.5 倍，或者 L1D 中热数据吞吐量的 3/4， 理论上，计数只是内部循环。这只是 inner-product 部分，未考虑不同的 pack/unpack 成本。
- 难以转置（但使用
pmovmskb to extract 1 bit from each byte and make them contiguous 可能并不可怕）。
每字节 6 个元素，0xxx0xxx（对于 HSW/SKL 上的这个问题可能没有优势，但值得考虑）
- 独立字节密度的 6 倍
- 通过 shifting/ORing 以交错方式从 0/1 字节创建相当容易，与每字节 4 位格式相同。
针对 AVX2 的高效 popcount 进行了优化 vpshufb。 2xvpshufb前无需屏蔽，1right-shift即可。 （如果设置了高位，vpshufb 将字节归零，否则它使用低半字节作为索引。这就是它需要屏蔽的原因。）将此格式右移 4 (vpsrld ymm0,4) 仍然会在每个字节的高位留零。 Load+AND -> accumulate popcount is 7 fused-domain uops per vector (vmovdqa/vpand ymm,[mem]/vpsrld ymm,4/2xvpshufb/2xvpaddb) ，其中只有 6 个需要 ALU 端口。因此 HSW/SKL 吞吐量在理论上是每 2 个时钟 1 个向量（192 个元素），或每个时钟 96 个元素。这需要每个时钟一个 256b 向量的平均负载吞吐量，因此正好可以解决 L2 带宽瓶颈。
理论上它与完全打包相同，但实际上它可能会稍微快一点或慢一点，具体取决于哪个安排更好（例如，更少的 AND/ADD uops 从洗牌中窃取端口 5）。完全打包可能更接近理论速度，因为它的更多 uops 可以 运行 在多个端口上。 Out-of-order 日程安排不完善的可能性较小。
pmovmskb 转置技巧不能正常工作。- 如果我们只需要
popcount(A[]) 而不是 popcount(A[] & B[])，这可能会有用。或者对于不同的微体系结构，其中 ALU 与负载吞吐量不同。
另一个变体，每个字节 7 个元素可以用单个 AVX512VBMI（Cannonlake？）vpermi2b (_mm512_permutex2var_epi8) 计算，其中每个索引字节选择 128 个字节之一来自其他两个寄存器的串联。如此宽的洗牌可能会很慢，但它有望比 AVX512 vpshufb separate-nibble 具有更好的吞吐量。
要使用 AVX512VBMI（但没有 AVX512VPOPCNTDQ）来计算 packed-8，您可以使用 vpermi2b 来计算低 7，然后 shift+mask 最高位并添加它。 （单个位的弹出计数 = 该位）。
uint8_t 元素更容易有效地洗牌（因为有像 vpshufb 这样的字节洗牌），所以如果你必须动态转置可能值得考虑。或者只在转置时即时打包成位？
32 位整数也是一种选择，但不是一个好的选择。每个向量的元素更少意味着转置中的洗牌指令更少，但不是原来的 4 倍。转置中的洗牌数量可能与 log2（每个向量的元素）之类的东西成比例。
这对于缓存占用空间/内存带宽来说也是一个大问题。 8 个大小差异的因子可能意味着做一整行或整列只占用 L1 的一部分，而不是溢出 L1。所以它可以使 cache-blocking 更容易/不那么重要。
10k * 20k / 8 = 每个矩阵 23.84MiB，使用 packed-bit 个元素。这比 L2 缓存大得多（在 Haswell 上为 256kiB，1MiB on Skylake-AVX512), but will fit in L3 on many-core Xeon CPUs. But L3 is competitively shared by all cores (including other VMs in a cloud environment), and is much slower than L2. (Many-core Xeons like you will be running on in HPC / cloud systems have lower per-core memory bandwidth than quad-core desktops, because of higher latency to L3 cache with no increase in concurrency (see the "latency-bound platforms" section of this answer。在 Xeon 上驱动相同数量的内存带宽需要更多的内核，即使总吞吐量更高。但是如果你能让每个内核大部分时间都在工作从它的私有 L2 中，你获得了很多。）
将 AND 结果相加：您已经安排了循环，因此您需要将单个 运行 布尔值减少到 non-zeros.这是好事。
对于 8 位整数 0/1 元素，在元素溢出之前最多可以执行 255 vpaddb。它具有良好的吞吐量：在 Haswell 上每个时钟 2 个，在 Skylake 上每个时钟 3 个。对于多个累加器，它涵盖了很多 AND 结果的向量。使用 vpsadbw against an all-zero vector to horizontally add the bytes in a vector into 64-bit integers. Then combine your accumulators with vpaddq, then horizontally sum it.
对于压缩位，您只想对 AND 结果的向量进行弹出计数。使用 AVX2 并且您的数据已经在向量中，您肯定想使用
VPSHUFB 基于 bit-slicing 人口统计。 （例如，参见 http://wm.ite.pl/articles/sse-popcount.html。如果必须手动对其进行矢量化，则需要使用内在函数而不是 asm 来编写它。）
您可以考虑将数据打包为每字节 4 位，在低半字节中。 这意味着 vpshufb 可以计算 AND 的每个字节中的位结果，不需要任何移位/掩蔽。在内部循环中，您将有 2 个负载，vpand、vpshufb、vpaddb。通过适当的展开，这应该跟上每时钟 2x 32B 的 L1D 负载带宽，并使所有三个向量执行端口（在 Haswell 或 Skylake 上）饱和。每 128 或 255 个向量或其他东西打破它，用 vpsadbw/vpaddq 累积你的累加器的字节。 （但是对于 cache-blocking，您可能无论如何都想经常中断并做不同的工作）。 所以 inner-most 循环应该 运行 每字节 4 个元素 * 每个向量 32B = 每个时钟周期 128 个元素， 如果你可以安排它读取数据在 L1D 缓存中很热。预计 Haswell/Skylake 上的 L2 缓存的带宽约为一半，L3 缓存的带宽会更糟。
有 uint8_t 个元素为 0 或 1，你也许可以使用一些整数 multiply-add 指令。它们的设计有点奇怪，用于与 FP FMA 不同的use-cases。它们添加水平成对的乘法结果，产生更宽的元素。 VPMADDUBSW 将 8 位元素加宽到 16 位元素，并且可以很好地处理 0 和 1。由于每个元素只能在 0..2 范围内，您仍然可以 horizontal-sum 和 vpsadbw。但是，如果您打算 vpsadbw，这对您没有任何好处 vpand。仅当您想使用 vpaddw 在向量累加器中使用 16 位元素而不是跳出循环以避免字节溢出时，它才有用。 vpmaddubsw doesn't seem useful here, becausevpsadbw` 是水平添加字节的更好方法。
使用SSE/AVX可以有效地将0/1整数转换为位图：对于32位整数元素，vpslld ymm0, 31到left-shift 相关位到每个元素的顶部，然后 vmovmskps eax, ymm0 得到每个 32 位元素的高字节的 8 位掩码。对于 8 位整数元素，vpslld ymm0, 7 / vpmovmskb eax, ymm0 做同样的事情，但对于每个字节，产生一个 32 位整数位图结果。 （只有每个字节的符号位很重要，所以没有只有 8 位粒度的移位指令很好。您不需要对携带到下一个元素的位做任何事情。）
这不是立即使用向量的好方法，因为您最终会在整数寄存器中得到结果。这不是即时生成和使用的好格式，但它是最紧凑的格式，因此如果您可以以这种格式 long-term 保存矩阵，它就很有意义。 （如果您在加载它们时会受到内存带宽的限制。）
将 32 位整数转换为 8 位整数：一种方法是使用 2x vpackssdw + vpacksswb。因为那些在 128b 通道内运行，所以您的元素最终将被重新排序。但只要每个 row/column 的顺序相同就可以了。如果您想获取不是以 32 元素的倍数开始的 row/column 块，这只是一个问题。这里的另一个选择是 left-shift（乘以 8、乘以 16 和乘以 24）和 OR 向量。实际上，您可以通过使用 1、2 或 3 个字节的未对齐加载偏移来免费进行移位。
static inline
__m256i load_interleave4x32(const int32_t *input) {
const char *p = (const char*)input;
__m256i t0 = _mm256_load_si256((const __m256i*)(p));
__m256i t1 = _mm256_load_si256((const __m256i*)(p+32*1-1)); // the 1/0 bits will be in the 2nd byte of each 32-bit element
__m256i t2 = _mm256_load_si256((const __m256i*)(p+32*2-2));
__m256i t3 = _mm256_load_si256((const __m256i*)(p+32*3-3));
return t0 | t1 | t2 | t3;
// or write this out with _mm256_or_si256, if you don't have overloaded operators like GNU C does.
// this should compile to 1 load and 3 vpor ymm0, [rdi+31] ... instructions.
}
转换为half-packed每字节4位：我们可以使用与上面相同的想法。从 load_interleave4x32 中获取 4 个向量（如果您从 8 位元素开始，则从 uint8_t 的数组中获取）。 Left-shift 它们按 0、1、2 和 3 位，然后将它们全部或。当我们只需要 AND a row/column 并弹出计算整个结果时，这个交错的 bit-order 很好，因为顺序无关紧要。 bit-order 解压缩回 in-order 字节相当有效，例如AND set1_epi8(1) 将为您提供一个字节向量。
如果您以这种格式存储整个矩阵，则可以将其用作转置的一部分，或者您可以使用这种格式存储 cache-blocked 转置的临时副本。 matmul 多次触及每个 row/column，因此第一次制作紧凑格式可能值得做额外的工作，因为这可以让您在后续传递中为每个向量做 4 倍的工作。
使用 AVX512BW (Skylake-AVX512)
我们真的很想用向量而不是标量整数来做 AND 和 popcnt，因为向量的宽度是 AVX2 的两倍，所以它们比标量 popcnt 更靠前。 （即使 Skylake-AVX512 在 运行 宁 512b 指令时关闭端口 1 上的矢量 ALU（但不是标量））。
@Harold points out an interesting identity 这让我们可以进行 2/3 的向量弹出计数，但代价是额外的整数运算。
popcnt(a) + popcnt(b) + popcnt(c)
= popcnt(a ^ b ^ c) + 2 * popcnt((a ^ b) & c | (a & b))
a ^ b ^ c 和 (a ^ b) & c | (a & b) 可以用一个 vpternlogd each (since each have 3 boolean inputs). The 2* is free if we use a separate pre-shifted vpshufb LUT vector. See also this implementation that uses 30x vpternlogd + 1 vector popcnt to handle 16 vectors of 512b 完成，最后进行一些清理（只在循环内执行 16*popcnt 计数；其他一切都是链接）。
这很可能值得计算 fully-packed 每字节元素 8 位，并且使该格式对 AVX512 更具吸引力，与 less-dense 针对 popcount 优化的格式相比没有那么多 shifting/masking.
vpternlogd 也可用作转置的 bit-blend 指令，如果 byte-granularity VPBLENDMB zmm{k1}, zmm, zmm 不是 fine-enough 粒度。
对于某些 CPU 上的 AVX2，这可能是值得的，也许可以避免每 4 或 5 个向量弹出计数中的 1 个而不是 3 个中的 1 个？或者如果它只是增加了总的执行端口压力并且没有任何特定的瓶颈，它可能根本没有帮助。它对标量 popcnt 指令很有用（可能在没有 AVX2 的 CPU 上），因为它们在 Intel CPU 的单个端口上确实存在瓶颈。
我们可以将 uint8_t 布尔元素转换为 non-interleaved 位图，比 AVX2 稍微更有效（甚至不需要移位），反之效率更高。 Test-into-mask 或 compare-into-mask 针对 set1_epi8(1) 的向量都可以完成这项工作，从 64 字节的输入中生成 64 位的掩码。或者从 32 位整数开始，一次生成 16 位掩码。您可以使用 kunpck 指令有效地连接这些位
_mm512_test_epi8_mask (vptestmb) 很有趣：AND 两个向量在一起，并产生 byte-elements 的 mask-register 结果是 true/false。但这并不是我们真正想要的：如果我们要打包我们的位，我们希望将其作为输入矩阵的 pre-processing 步骤，而不是在进行内积时即时进行。
bitmap -> vector of 0 / -1 很快：__m512i _mm512_movm_epi8 (__mmask64 k) (vpmovm2b) 在一条指令中完成。您可以减去 -1 而不是添加 1，但您必须先屏蔽它，然后才能将一个字节内的多个位或运算在一起。
没有 AVX512BW 或 AVX512DQ（Knight's Landing Xeon Phi），您没有 512b vpshufb，因此您无法有效地向量 popcnt。有一个AVX512 popcnt extension for vector popcnt directly, but no hardware with it has even been announced yet. (AVX2 vpshufb ymm is very slow on KNL, though: one per 12 cycles, and psadbw ymm is 1 per 9 cycles, so even using 256b vectors is unattractive). You might use a bithack popcnt based on 32-bit integer elements, since that's just AND/shift/ADD。 32 位元素比 64 位元素需要更少的步骤来 popcnt，但是对于合理的问题大小来说仍然足够大而不会溢出（所以你可以将向量的水平和推迟到循环之外）
考虑到存储格式的选择，每个字节打包多个位对于 KNL 可能不是一个好主意，但 single-byte 整数元素是个好主意。 vpandd zmm 和 vpaddd zmm 都很快并且是 AVX512F 的一部分，我们可以使用它们，因为我们不想让我们的 single-byte 溢出。 （当我们实际上有 8 位元素不会相互携带时使用打包的 32 位加法是一种 SWAR 技术。）相对于 Skylake-AVX512，KNL 具有良好的内存带宽和较差的指令吞吐量，我想想。
转位：
BMI2 _pdep_u64 在这里可能会有用。这是一个标量 instruction/intrinsic。如果它使 bit-transpose 比解包到字节更有效，您可能希望在使用 AND + 计数的向量加载重新加载它之前存储一个转置结果块。 （在标量存储后立即重新加载向量将导致 store-forwarding 停顿。）
另一个有用的选项是 vpmovmskb 可以从 32 字节向量中切出 32 位，每个字节一个。这为您提供了转置的构建块，可能与字节混洗相结合以按正确的顺序获取字节。有关详细信息，请参阅 this blog post, and also How would you transpose a binary matrix?。
在 matmul 中使用它
您的某些选择取决于输入数据的格式，以及重复使用相同矩阵的频率。如果一个矩阵将被多次使用，那么提前将其压缩为每字节 4 或 8 位是有意义的。 （或者在第一次使用时即时）。保留它的转置副本也可能有意义，特别是如果它始终是需要转置的乘法的一侧。 （如果您有时需要一种方式，有时需要另一种方式，则动态重做可能更适合 L3 缓存占用空间。但是这些足够大，您可能不会获得很多 L3 命中，因此只需保留转置副本即可好。）
或者甚至可以在从您的输入格式转换时写出转置和 non-transposed 版本。
你肯定会想要 cache-block 相乘，所以相同的数据在 L1 中很热时会被多次重复使用。关于这件事，我没有什么要说的。 与 cache-blocking 普通 FP matmul 时相同的原则适用，所以去读一读吧。
对您的 C++ 实现的评论：
对整个列使用位集 & 会将值放回内存中，然后您将在结果的 .count() 中再次循环它们。我怀疑编译器会将其优化为 one-pass 循环，该循环在 VPAND 结果的每个向量上使用基于 VPSHUFB 的 bit-slicing popcnt，但这会好得多。 （例如，参见 http://wm.ite.pl/articles/sse-popcount.html。如果必须手动对其进行矢量化，则需要使用内在函数而不是 asm 来编写它。）
对于您的矩阵大小，至少内部循环可能会命中 L1D 缓存，但是循环两次的额外 load/store 指令会产生更多开销，并且还会干扰有价值数据的预取。
让编译器高效地弹出 dynamically-sized 位图（无需手动矢量化） 并不容易。唯一不烂的是 clang++ -stdlib=libc++ 和 vector<bool>，它将 std::count(v.begin(), v.end(), true); 编译为矢量化 vpshufb + vpsadbw + vpaddq 循环，这很好。如果它只在展开的循环中使用 vpaddb 并且每次迭代使用一次 vpsadbw + vpaddq 会更快，但它对于 auto-vectorized 代码非常好。
g++ 的 vector<bool> 也是一个位图，但是 std::count(v.begin(), v.end(), true); 非常糟糕：它使用了一个完全天真的循环，一次测试 1 位。它甚至不能有效地做到这一点。与 clang++ 相同，默认 libstdc++ 而不是新的 libc++.
boost::dynamic_bitset 有一个 .count() 成员函数，但它没有利用 popcnt 指令或 AVX2。它执行 byte-at-a-time LUT 查找。这比没有 libc++ 的 std::count(vector<bool>) 要好得多，但对于 HPC 来说还不够好。
这是测试代码 on the Godbolt compiler explorer，带有 gcc 和 clang asm 输出。都用了-march=haswell.
但不幸的是，似乎没有一种有效的方法来 bitwise-AND 两个std::vector<bool>。 This answer 展示了如何获得 g++ libstdc++ vector<bool> 的底层实现，但该代码没有 auto-vectorize。对 libc++ 做同样的事情并对其进行调整，使其 auto-vectorizes 可能 让您获得手动矢量化可能的大部分性能（转置除外） ，但您可能必须将整个矩阵保持在一个 vector<bool> 中，因为向量的向量是一个糟糕的额外间接级别。如果问题的转置部分也是 performance-critical，使用标准容器来访问有效的 popcount 并不能解决整个问题。
对于 std::bitset<1024*1024>.count()，clang 在有或没有 libc++ 的情况下都能产生同样有效的 AVX2 popcount。 g++ 使用 64 位 popcnt 指令进行标量循环，在 Haswell 和 Skylake 上，该指令（根据 this）对于小比特集比良好的 AVX2 popcnt 要快一些，但对于大比特集来说要慢一些.
另请参阅：On vector<bool> — Howard Hinnant，了解有关 C++ 标准库的一些评论，以及为什么位数组是有用的数据结构，但 vector<bool> 是一个不好的名称。此外，properly-optimized count/find_first/etc 的一些基准测试。在 bit-vector 与 1 bool-per-byte bool[] 数组上，与天真的 vector<bool> 相比（就像你从没有 libc++ 的 gcc 和 clang 得到的一样）。
对于大型二进制矩阵（10Kx20K）的乘法，我通常做的是将矩阵转换为浮点矩阵并执行浮点矩阵乘法，因为整数矩阵乘法非常慢（
不过这一次，我需要执行数十万次这样的乘法运算，平均而言，即使是毫秒级的性能提升对我来说也很重要。
我想要一个 int 或 float 矩阵作为结果，因为乘积可能包含非 0 或 1 的元素。输入矩阵元素全部为 0 或 1，因此它们可以存储为单个位。
在行向量和列向量之间的内积（产生输出矩阵的一个元素）中，乘法简化为按位与。加法仍然是加法，但我们可以使用人口计数函数添加位，而不是单独循环它们。
一些其他 boolean/binary-matrix 函数或位而不是计算它们，产生位矩阵结果，但这不是我需要的。
这是一个示例代码，显示将问题形成为 std::bitset、AND 和 count 运算比矩阵乘法更快。
#include <iostream>
using std::cout; using std::endl;
#include <vector>
using std::vector;
#include <chrono>
#include <Eigen/Dense>
using Eigen::Map; using Eigen::Matrix; using Eigen::MatrixXf;
#include <random>
using std::random_device; using std::mt19937; using std::uniform_int_distribution;
#include <bitset>
using std::bitset;
using std::floor;
const int NROW = 1000;
const int NCOL = 20000;
const float DENSITY = 0.4;
const float DENOMINATOR = 10.0 - (10*DENSITY);
void fill_random(vector<float>& vec) {
random_device rd;
mt19937 eng(rd());
uniform_int_distribution<> distr(0, 10);
int nnz = 0;
for (int i = 0; i < NROW*NCOL; ++i)
vec.push_back(floor(distr(eng)/DENOMINATOR));
}
void matmul(vector<float>& vec){
float *p = vec.data();
MatrixXf A = Eigen::Map<Eigen::Matrix<float, NROW, NCOL, Eigen::RowMajor>>(p);
cout << "Eigen matrix has " << A.rows() << " rows and " << A.cols() << " columns." << endl;
cout << "Total non-zero values : " << A.sum() << endl;
cout << "The density of non-zero values is " << A.sum() * 1.0 / (A.cols()*A.rows()) << endl;
auto start = std::chrono::steady_clock::now();
MatrixXf B = A.transpose() * A;
auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
cout << "Mat mul took " << end << " ms"<< endl;
// Just to make sure the operation is not skipped by compiler
cout << "Eigen coo ";
for (int i=0; i<10; ++i)
cout << B(0,i) << " ";
cout << endl;
}
void bitset_op(vector<float>& vec) {
// yeah it's not a great idea to set size at compile time but have to
vector<bitset<NROW>> col_major(NCOL);
// right, multiple par for isn't a good idea, maybe in a parallel block
// Doing this for simplicity to profile second loop timing
// converting row major float vec to col major bool vec
#pragma omp parallel for
for (int j=0; j < NCOL; ++j) {
for (int i=0; i < NROW; ++i) {
col_major[j].set(i, vec[i*NCOL + j] && 1);
}
}
auto start = std::chrono::steady_clock::now();
vector<int> coo;
coo.assign(NCOL*NCOL, 0);
#pragma omp parallel for
for (int j=0; j < NCOL; ++j) {
for (int k=0; k<NCOL; ++k) {
coo[j*NCOL + k] = (col_major[j]&col_major[k]).count();
}
}
auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
cout << "bitset intersection took " << end << " ms"<< endl;
// Just to make sure the operation is not skipped by compiler
cout << "biset coo ";
for (int i=0; i<10; ++i)
cout << coo[i] << " ";
cout << endl;
}
int main() {
// Saving to float instead of int to speed up matmul
vector<float> vec;
fill_random(vec);
matmul(vec);
bitset_op(vec);
}
运行 这与：
g++ -O3 -fopenmp -march=native -I. -std=c++11 code.cpp -o code
我得到：
Eigen matrix has 1000 rows and 20000 columns.
Total non-zero values : 9.08978e+06
The density of non-zero values is 0.454489
Mat mul took 1849 ms
Eigen coo 458 206 208 201 224 205 204 199 217 210
bitset intersection took 602 ms
biset coo 458 206 208 201 224 205 204 199 217 210
如您所见，matmul 作为一组 bitset 操作比 Eigen 的 float matmul 快大约 3 倍，这是有道理的。
我想强调的是，我需要在 100K（在 HPC 或云中）上执行此操作，平均毫秒级的性能提升会有所不同。
我不受任何特定库、C++ 标准等的约束。所以请随时回答您认为比使用 GPU 的解决方案更快的任何解决方案，因为我不能将它用于数字原因。
我不确定你可以让编译器为你做多少，如果有的话，而无需使用内在函数或 C++ vector-class 包装器手动矢量化（如 Agner Fog's VCL，如果你的项目是许可证与 GPL 兼容）。还有一些 non-GPLed 包装器。
Cache-blocking 矩阵乘法是一门艺术（在这里很重要），如果你可以使用 Eigen 的现有模板，那将是非常好的，但是一个不同的 class，它在整数上使用按位 and，而不是在浮点数上乘法。我不确定这是否可行。
我做了一些搜索，大多数关于二进制矩阵的文献都是关于产生布尔结果的（包括 SO 问题 like this）。向量 inner-product 是用 AND 作为乘法，但用 XOR 或 OR 作为加法，而不是 popcount。也许有一个 search-term 我错过了它描述 "normal" 矩阵恰好是 (0,1) 矩阵，但乘积不是。
由于每一毫秒都很重要，您可能需要手动对其进行矢量化。
并不是说 vector-integer 东西一般都很慢，只是 只是 vector-integer 相乘比较慢，尤其是与 vector-float 相比最近的 x86 硬件上的 FMA（尤其是英特尔，它在 Haswell 及更高版本上的 FP FMA 吞吐量为每时钟 2x 256b 矢量）。
由于您不需要实际与布尔元素相乘，只需一个 AND（每个时钟吞吐量 3 个向量），这对您来说不是问题。 efficiency-gain 从每个向量中获取更多元素应该足以弥补任何额外成本 per-vector.
当然，这假设一个整数 matmul 实现使用所有相同的 cache-blocking 和其他优化作为等效的 FP matmul，如果你不想（或不想），这就是问题所在知道如何）自己编写，找不到可以为您编写的库。
我只是在回答一个问题，即通过最佳实施， 的效率如何。 标题问题的答案是非常明确的是，使用实际乘法是浪费大量时间，尤其是对于 32 位元素。
存储格式选项：
每字节一个元素(0/1):
- 4 倍密度
float（缓存占用空间/内存带宽/每个向量的元素） - 易于使用字节洗牌进行转置
- 垂直 ADD 很容易，以防万一（例如，在外循环上进行矢量化，并同时处理多行或多列。如果您有数据，可能会很好（避免最后水平求和）交错的方式使这项工作无需额外的改组。）
每字节4个元素，打包到低半字节:
- 4 倍的单独字节密度
- 使用 AVX2
vpshufb进行 popcount 非常有效。对于 L1D 缓存中的热输入，理论上，您可以 load/AND/accumulate-a-popcount 每个时钟周期（每个内核）吞吐 128 AND-result 个元素。每个时钟 4 fused-domain 微指令饱和 SKL/HSW front-end 发出每个时钟 4 的带宽，并且不会在 3 个矢量 ALU 端口上出现瓶颈，因为其中一个微指令是纯负载。 （另一个负载 micro-fuses 与vpand）。当 L2 带宽出现瓶颈时（每个周期约 32B 负载），运行s 在每个时钟 64 个元素。见下文。 - 从整数或 packed-bitmap 创建速度较慢（但如果您以交错顺序将位放入向量以实现高效 pack/unpack 到 in-order 字节，而不是强制位到井然有序）。
- 难以转置（可能比满载还差）
压缩位:
- 8 倍单独字节的密度，每个 AVX2 向量 256 个元素。
- 可以从具有
pmovmskb的向量创建，用于 non-interleaved 存储顺序。 （不过，对于动态创建不是很有用，因为这会将结果放在整数 reg 中，而不是向量中。交错的 bit-order 可能是最好的，尤其是在转置期间解包）。 - 使用 AVX2 进行 popcount 相当有效：mask / shift+mask / 2x
vpshufb。 (9 fused-domain uops (8 vector-ALU uops) 到 AND + 累积 256 个元素的 popcount（来自 2 row/column 向量），对比 8 uops (6 vector-ALU uops) 4-per-byte 策略（来自 4 row/column 向量）。ALU 端口瓶颈将其限制为来自 L1D 或 L2 的每个时钟 96 个元素。因此，当它在 L2 带宽上出现瓶颈时，这大约是 pack4 策略的 inner-product 吞吐量的 1.5 倍，或者 L1D 中热数据吞吐量的 3/4， 理论上，计数只是内部循环。这只是 inner-product 部分，未考虑不同的 pack/unpack 成本。 - 难以转置（但使用
pmovmskbto extract 1 bit from each byte and make them contiguous 可能并不可怕）。
每字节 6 个元素，0xxx0xxx（对于 HSW/SKL 上的这个问题可能没有优势，但值得考虑）
- 独立字节密度的 6 倍
- 通过 shifting/ORing 以交错方式从 0/1 字节创建相当容易，与每字节 4 位格式相同。
针对 AVX2 的高效 popcount 进行了优化
vpshufb。 2xvpshufb前无需屏蔽，1right-shift即可。 （如果设置了高位，vpshufb将字节归零，否则它使用低半字节作为索引。这就是它需要屏蔽的原因。）将此格式右移 4 (vpsrld ymm0,4) 仍然会在每个字节的高位留零。 Load+AND -> accumulate popcount is 7 fused-domain uops per vector (vmovdqa/vpand ymm,[mem]/vpsrld ymm,4/2xvpshufb/2xvpaddb) ，其中只有 6 个需要 ALU 端口。因此 HSW/SKL 吞吐量在理论上是每 2 个时钟 1 个向量（192 个元素），或每个时钟 96 个元素。这需要每个时钟一个 256b 向量的平均负载吞吐量，因此正好可以解决 L2 带宽瓶颈。理论上它与完全打包相同，但实际上它可能会稍微快一点或慢一点，具体取决于哪个安排更好（例如，更少的 AND/ADD uops 从洗牌中窃取端口 5）。完全打包可能更接近理论速度，因为它的更多 uops 可以 运行 在多个端口上。 Out-of-order 日程安排不完善的可能性较小。
pmovmskb转置技巧不能正常工作。- 如果我们只需要
popcount(A[])而不是popcount(A[] & B[])，这可能会有用。或者对于不同的微体系结构，其中 ALU 与负载吞吐量不同。
另一个变体，每个字节 7 个元素可以用单个 AVX512VBMI（Cannonlake？）vpermi2b (_mm512_permutex2var_epi8) 计算，其中每个索引字节选择 128 个字节之一来自其他两个寄存器的串联。如此宽的洗牌可能会很慢，但它有望比 AVX512 vpshufb separate-nibble 具有更好的吞吐量。
要使用 AVX512VBMI（但没有 AVX512VPOPCNTDQ）来计算 packed-8，您可以使用 vpermi2b 来计算低 7，然后 shift+mask 最高位并添加它。 （单个位的弹出计数 = 该位）。
uint8_t 元素更容易有效地洗牌（因为有像 vpshufb 这样的字节洗牌），所以如果你必须动态转置可能值得考虑。或者只在转置时即时打包成位？
32 位整数也是一种选择，但不是一个好的选择。每个向量的元素更少意味着转置中的洗牌指令更少，但不是原来的 4 倍。转置中的洗牌数量可能与 log2（每个向量的元素）之类的东西成比例。
这对于缓存占用空间/内存带宽来说也是一个大问题。 8 个大小差异的因子可能意味着做一整行或整列只占用 L1 的一部分，而不是溢出 L1。所以它可以使 cache-blocking 更容易/不那么重要。
10k * 20k / 8 = 每个矩阵 23.84MiB，使用 packed-bit 个元素。这比 L2 缓存大得多（在 Haswell 上为 256kiB，1MiB on Skylake-AVX512), but will fit in L3 on many-core Xeon CPUs. But L3 is competitively shared by all cores (including other VMs in a cloud environment), and is much slower than L2. (Many-core Xeons like you will be running on in HPC / cloud systems have lower per-core memory bandwidth than quad-core desktops, because of higher latency to L3 cache with no increase in concurrency (see the "latency-bound platforms" section of this answer。在 Xeon 上驱动相同数量的内存带宽需要更多的内核，即使总吞吐量更高。但是如果你能让每个内核大部分时间都在工作从它的私有 L2 中，你获得了很多。）
将 AND 结果相加：您已经安排了循环，因此您需要将单个 运行 布尔值减少到 non-zeros.这是好事。
对于 8 位整数 0/1 元素，在元素溢出之前最多可以执行 255 vpaddb。它具有良好的吞吐量：在 Haswell 上每个时钟 2 个，在 Skylake 上每个时钟 3 个。对于多个累加器，它涵盖了很多 AND 结果的向量。使用 vpsadbw against an all-zero vector to horizontally add the bytes in a vector into 64-bit integers. Then combine your accumulators with vpaddq, then horizontally sum it.
对于压缩位，您只想对 AND 结果的向量进行弹出计数。使用 AVX2 并且您的数据已经在向量中，您肯定想使用
VPSHUFB 基于 bit-slicing 人口统计。 （例如，参见 http://wm.ite.pl/articles/sse-popcount.html。如果必须手动对其进行矢量化，则需要使用内在函数而不是 asm 来编写它。）
您可以考虑将数据打包为每字节 4 位，在低半字节中。 这意味着 vpshufb 可以计算 AND 的每个字节中的位结果，不需要任何移位/掩蔽。在内部循环中，您将有 2 个负载，vpand、vpshufb、vpaddb。通过适当的展开，这应该跟上每时钟 2x 32B 的 L1D 负载带宽，并使所有三个向量执行端口（在 Haswell 或 Skylake 上）饱和。每 128 或 255 个向量或其他东西打破它，用 vpsadbw/vpaddq 累积你的累加器的字节。 （但是对于 cache-blocking，您可能无论如何都想经常中断并做不同的工作）。 所以 inner-most 循环应该 运行 每字节 4 个元素 * 每个向量 32B = 每个时钟周期 128 个元素， 如果你可以安排它读取数据在 L1D 缓存中很热。预计 Haswell/Skylake 上的 L2 缓存的带宽约为一半，L3 缓存的带宽会更糟。
有 uint8_t 个元素为 0 或 1，你也许可以使用一些整数 multiply-add 指令。它们的设计有点奇怪，用于与 FP FMA 不同的use-cases。它们添加水平成对的乘法结果，产生更宽的元素。 VPMADDUBSW 将 8 位元素加宽到 16 位元素，并且可以很好地处理 0 和 1。由于每个元素只能在 0..2 范围内，您仍然可以 horizontal-sum 和 vpsadbw。但是，如果您打算 vpsadbw，这对您没有任何好处 vpand。仅当您想使用 vpaddw 在向量累加器中使用 16 位元素而不是跳出循环以避免字节溢出时，它才有用。 vpmaddubsw doesn't seem useful here, becausevpsadbw` 是水平添加字节的更好方法。
使用SSE/AVX可以有效地将0/1整数转换为位图：对于32位整数元素，vpslld ymm0, 31到left-shift 相关位到每个元素的顶部，然后 vmovmskps eax, ymm0 得到每个 32 位元素的高字节的 8 位掩码。对于 8 位整数元素，vpslld ymm0, 7 / vpmovmskb eax, ymm0 做同样的事情，但对于每个字节，产生一个 32 位整数位图结果。 （只有每个字节的符号位很重要，所以没有只有 8 位粒度的移位指令很好。您不需要对携带到下一个元素的位做任何事情。）
这不是立即使用向量的好方法，因为您最终会在整数寄存器中得到结果。这不是即时生成和使用的好格式，但它是最紧凑的格式，因此如果您可以以这种格式 long-term 保存矩阵，它就很有意义。 （如果您在加载它们时会受到内存带宽的限制。）
将 32 位整数转换为 8 位整数：一种方法是使用 2x vpackssdw + vpacksswb。因为那些在 128b 通道内运行，所以您的元素最终将被重新排序。但只要每个 row/column 的顺序相同就可以了。如果您想获取不是以 32 元素的倍数开始的 row/column 块，这只是一个问题。这里的另一个选择是 left-shift（乘以 8、乘以 16 和乘以 24）和 OR 向量。实际上，您可以通过使用 1、2 或 3 个字节的未对齐加载偏移来免费进行移位。
static inline
__m256i load_interleave4x32(const int32_t *input) {
const char *p = (const char*)input;
__m256i t0 = _mm256_load_si256((const __m256i*)(p));
__m256i t1 = _mm256_load_si256((const __m256i*)(p+32*1-1)); // the 1/0 bits will be in the 2nd byte of each 32-bit element
__m256i t2 = _mm256_load_si256((const __m256i*)(p+32*2-2));
__m256i t3 = _mm256_load_si256((const __m256i*)(p+32*3-3));
return t0 | t1 | t2 | t3;
// or write this out with _mm256_or_si256, if you don't have overloaded operators like GNU C does.
// this should compile to 1 load and 3 vpor ymm0, [rdi+31] ... instructions.
}
转换为half-packed每字节4位：我们可以使用与上面相同的想法。从 load_interleave4x32 中获取 4 个向量（如果您从 8 位元素开始，则从 uint8_t 的数组中获取）。 Left-shift 它们按 0、1、2 和 3 位，然后将它们全部或。当我们只需要 AND a row/column 并弹出计算整个结果时，这个交错的 bit-order 很好，因为顺序无关紧要。 bit-order 解压缩回 in-order 字节相当有效，例如AND set1_epi8(1) 将为您提供一个字节向量。
如果您以这种格式存储整个矩阵，则可以将其用作转置的一部分，或者您可以使用这种格式存储 cache-blocked 转置的临时副本。 matmul 多次触及每个 row/column，因此第一次制作紧凑格式可能值得做额外的工作，因为这可以让您在后续传递中为每个向量做 4 倍的工作。
使用 AVX512BW (Skylake-AVX512)
我们真的很想用向量而不是标量整数来做 AND 和 popcnt，因为向量的宽度是 AVX2 的两倍，所以它们比标量 popcnt 更靠前。 （即使 Skylake-AVX512 在 运行 宁 512b 指令时关闭端口 1 上的矢量 ALU（但不是标量））。
@Harold points out an interesting identity 这让我们可以进行 2/3 的向量弹出计数，但代价是额外的整数运算。
popcnt(a) + popcnt(b) + popcnt(c)
= popcnt(a ^ b ^ c) + 2 * popcnt((a ^ b) & c | (a & b))
a ^ b ^ c 和 (a ^ b) & c | (a & b) 可以用一个 vpternlogd each (since each have 3 boolean inputs). The 2* is free if we use a separate pre-shifted vpshufb LUT vector. See also this implementation that uses 30x vpternlogd + 1 vector popcnt to handle 16 vectors of 512b 完成，最后进行一些清理（只在循环内执行 16*popcnt 计数；其他一切都是链接）。
这很可能值得计算 fully-packed 每字节元素 8 位，并且使该格式对 AVX512 更具吸引力，与 less-dense 针对 popcount 优化的格式相比没有那么多 shifting/masking.
vpternlogd 也可用作转置的 bit-blend 指令，如果 byte-granularity VPBLENDMB zmm{k1}, zmm, zmm 不是 fine-enough 粒度。
对于某些 CPU 上的 AVX2，这可能是值得的，也许可以避免每 4 或 5 个向量弹出计数中的 1 个而不是 3 个中的 1 个？或者如果它只是增加了总的执行端口压力并且没有任何特定的瓶颈，它可能根本没有帮助。它对标量 popcnt 指令很有用（可能在没有 AVX2 的 CPU 上），因为它们在 Intel CPU 的单个端口上确实存在瓶颈。
我们可以将 uint8_t 布尔元素转换为 non-interleaved 位图，比 AVX2 稍微更有效（甚至不需要移位），反之效率更高。 Test-into-mask 或 compare-into-mask 针对 set1_epi8(1) 的向量都可以完成这项工作，从 64 字节的输入中生成 64 位的掩码。或者从 32 位整数开始，一次生成 16 位掩码。您可以使用 kunpck 指令有效地连接这些位
_mm512_test_epi8_mask (vptestmb) 很有趣：AND 两个向量在一起，并产生 byte-elements 的 mask-register 结果是 true/false。但这并不是我们真正想要的：如果我们要打包我们的位，我们希望将其作为输入矩阵的 pre-processing 步骤，而不是在进行内积时即时进行。
bitmap -> vector of 0 / -1 很快：__m512i _mm512_movm_epi8 (__mmask64 k) (vpmovm2b) 在一条指令中完成。您可以减去 -1 而不是添加 1，但您必须先屏蔽它，然后才能将一个字节内的多个位或运算在一起。
没有 AVX512BW 或 AVX512DQ（Knight's Landing Xeon Phi），您没有 512b vpshufb，因此您无法有效地向量 popcnt。有一个AVX512 popcnt extension for vector popcnt directly, but no hardware with it has even been announced yet. (AVX2 vpshufb ymm is very slow on KNL, though: one per 12 cycles, and psadbw ymm is 1 per 9 cycles, so even using 256b vectors is unattractive). You might use a bithack popcnt based on 32-bit integer elements, since that's just AND/shift/ADD。 32 位元素比 64 位元素需要更少的步骤来 popcnt，但是对于合理的问题大小来说仍然足够大而不会溢出（所以你可以将向量的水平和推迟到循环之外）
考虑到存储格式的选择，每个字节打包多个位对于 KNL 可能不是一个好主意，但 single-byte 整数元素是个好主意。 vpandd zmm 和 vpaddd zmm 都很快并且是 AVX512F 的一部分，我们可以使用它们，因为我们不想让我们的 single-byte 溢出。 （当我们实际上有 8 位元素不会相互携带时使用打包的 32 位加法是一种 SWAR 技术。）相对于 Skylake-AVX512，KNL 具有良好的内存带宽和较差的指令吞吐量，我想想。
转位：
BMI2 _pdep_u64 在这里可能会有用。这是一个标量 instruction/intrinsic。如果它使 bit-transpose 比解包到字节更有效，您可能希望在使用 AND + 计数的向量加载重新加载它之前存储一个转置结果块。 （在标量存储后立即重新加载向量将导致 store-forwarding 停顿。）
另一个有用的选项是 vpmovmskb 可以从 32 字节向量中切出 32 位，每个字节一个。这为您提供了转置的构建块，可能与字节混洗相结合以按正确的顺序获取字节。有关详细信息，请参阅 this blog post, and also How would you transpose a binary matrix?。
在 matmul 中使用它
您的某些选择取决于输入数据的格式，以及重复使用相同矩阵的频率。如果一个矩阵将被多次使用，那么提前将其压缩为每字节 4 或 8 位是有意义的。 （或者在第一次使用时即时）。保留它的转置副本也可能有意义，特别是如果它始终是需要转置的乘法的一侧。 （如果您有时需要一种方式，有时需要另一种方式，则动态重做可能更适合 L3 缓存占用空间。但是这些足够大，您可能不会获得很多 L3 命中，因此只需保留转置副本即可好。）
或者甚至可以在从您的输入格式转换时写出转置和 non-transposed 版本。
你肯定会想要 cache-block 相乘，所以相同的数据在 L1 中很热时会被多次重复使用。关于这件事，我没有什么要说的。 与 cache-blocking 普通 FP matmul 时相同的原则适用，所以去读一读吧。
对您的 C++ 实现的评论：
对整个列使用位集 & 会将值放回内存中，然后您将在结果的 .count() 中再次循环它们。我怀疑编译器会将其优化为 one-pass 循环，该循环在 VPAND 结果的每个向量上使用基于 VPSHUFB 的 bit-slicing popcnt，但这会好得多。 （例如，参见 http://wm.ite.pl/articles/sse-popcount.html。如果必须手动对其进行矢量化，则需要使用内在函数而不是 asm 来编写它。）
对于您的矩阵大小，至少内部循环可能会命中 L1D 缓存，但是循环两次的额外 load/store 指令会产生更多开销，并且还会干扰有价值数据的预取。
让编译器高效地弹出 dynamically-sized 位图（无需手动矢量化） 并不容易。唯一不烂的是 clang++ -stdlib=libc++ 和 vector<bool>，它将 std::count(v.begin(), v.end(), true); 编译为矢量化 vpshufb + vpsadbw + vpaddq 循环，这很好。如果它只在展开的循环中使用 vpaddb 并且每次迭代使用一次 vpsadbw + vpaddq 会更快，但它对于 auto-vectorized 代码非常好。
g++ 的 vector<bool> 也是一个位图，但是 std::count(v.begin(), v.end(), true); 非常糟糕：它使用了一个完全天真的循环，一次测试 1 位。它甚至不能有效地做到这一点。与 clang++ 相同，默认 libstdc++ 而不是新的 libc++.
boost::dynamic_bitset 有一个 .count() 成员函数，但它没有利用 popcnt 指令或 AVX2。它执行 byte-at-a-time LUT 查找。这比没有 libc++ 的 std::count(vector<bool>) 要好得多，但对于 HPC 来说还不够好。
这是测试代码 on the Godbolt compiler explorer，带有 gcc 和 clang asm 输出。都用了-march=haswell.
但不幸的是，似乎没有一种有效的方法来 bitwise-AND 两个std::vector<bool>。 This answer 展示了如何获得 g++ libstdc++ vector<bool> 的底层实现，但该代码没有 auto-vectorize。对 libc++ 做同样的事情并对其进行调整，使其 auto-vectorizes 可能 让您获得手动矢量化可能的大部分性能（转置除外） ，但您可能必须将整个矩阵保持在一个 vector<bool> 中，因为向量的向量是一个糟糕的额外间接级别。如果问题的转置部分也是 performance-critical，使用标准容器来访问有效的 popcount 并不能解决整个问题。
对于 std::bitset<1024*1024>.count()，clang 在有或没有 libc++ 的情况下都能产生同样有效的 AVX2 popcount。 g++ 使用 64 位 popcnt 指令进行标量循环，在 Haswell 和 Skylake 上，该指令（根据 this）对于小比特集比良好的 AVX2 popcnt 要快一些，但对于大比特集来说要慢一些.
另请参阅：On vector<bool> — Howard Hinnant，了解有关 C++ 标准库的一些评论，以及为什么位数组是有用的数据结构，但 vector<bool> 是一个不好的名称。此外，properly-optimized count/find_first/etc 的一些基准测试。在 bit-vector 与 1 bool-per-byte bool[] 数组上，与天真的 vector<bool> 相比（就像你从没有 libc++ 的 gcc 和 clang 得到的一样）。