代码之家 › 专栏 › 技术社区 › NULL

使用按位AND和popcount而不是实际整数或浮点乘法的大(0,1)矩阵乘法?

bitset avx matrix-multiplication sse c++

10

NULL · 技术社区 · 7 年前

have a look at here

不过,这次我需要执行数十万次乘法 对我来说,平均而言,即使是毫秒级的性能提升也很重要 .

int float 结果是矩阵,因为乘积可能有非0或1的元素。输入矩阵元素都是0或1,因此它们可以存储为单个位。

在行向量和列向量之间的内积中(生成输出矩阵的一个元素),乘法简化为按位and。加法仍然是加法,但我们可以使用人口计数函数来添加位,而不是逐个循环。

其他一些布尔/二进制矩阵函数或位,而不是计算它们,生成位矩阵结果,但这不是我需要的。

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大约比Eigen的浮点matmul快3倍,这是合理的。

我想强调的是,我需要在10万以上执行此操作

我不受任何特定库、C++标准等的约束。因此,请随意回答任何您认为比使用GPU更快的解决方案,因为我有很多原因不能使用它。

1 回复 | 直到 7 年前

1

16

Peter Cordes 7 年前

Agner Fog's VCL ,如果您的项目许可证与GPL兼容)。还有一些非GPLed包装。

,如果您可以使用Eigen的现有模板,但使用另一个使用位的类,那将非常好 and 在整数上,而不是在浮点上乘法。我不确定这是否可行。

我做了一些搜索,大多数关于二进制矩阵的文献都是关于生成布尔结果的(包括一些问题) like this

一般来说,向量整数并不是很慢,而是只是向量整数乘法速度很慢,尤其是与向量相比- float 最新x86硬件上的FMA(特别是Intel,其在Haswell和更高版本上的浮点FMA吞吐量为每个时钟2x 256b向量)。

当然,这假设一个整数matmul实现使用与等效FP matmul相同的缓存阻塞和其他优化,如果您不想(或不知道如何)自己编写,并且找不到一个库来为您编写,那么问题就出在这里。

标题问题的答案是非常明确的 ,使用实际乘法非常浪费时间,尤其是使用32位元素。

存储格式选项:

密度的4倍 (每个向量的缓存占用空间/内存带宽/元素)
在重要的情况下,垂直加法很容易(例如,在外循环上进行矢量化,并同时处理多行或多列。如果您以一种不需要额外洗牌的方式交错数据,则可以很好(避免最后的水平和)

每字节4个元素,压缩到低半字节中 :

使用AVX2进行popcount非常有效 vpshufb vpand ). 当L2带宽受限时(每个周期约32B负载),每个时钟运行64个元素。见下文。

压缩位 :

,每个AVX2矢量256个元素。
可以从矢量创建 pmovmskb
使用AVX2对popcount相当有效:掩码/移位+掩码/2x vpshufb (9个融合域UOP(8个向量ALU UOP)和+累积256个元素的popcount(来自2行/列向量),而8个UOP(6个向量ALU UOP)用于每字节4个策略(来自4行/列向量)ALU端口瓶颈将L1D或L2的每个时钟限制为96个元素 当二级带宽出现瓶颈时,这大约是pack4策略的内积吞吐量的1.5倍,或者是L1D中数据热的吞吐量的3/4,
很难转换(但可能不可怕 pmovmskb to extract 1 bit from each byte and make them contiguous

每字节6个元素, 0xxx0xxx (这个问题在HSW/SKL上可能没有优势,但值得考虑):

6x独立字节的密度
使用AVX2优化了有效的popcount vpshufb 2倍前无需遮罩 ,只需右移一次。( 如果设置了高位,则将字节归零,否则将使用低位半字节作为索引。这就是为什么它需要屏蔽。)将此格式右移4( vpsrld ymm0,4 )仍将在每个字节的高位保留零。加载+和->累积popcount是每个向量7个融合域UOP( vmovdqa vpand ymm,[mem] / vpsrld ymm,4 vpshufb /2倍 vpaddb
这个
如果我们需要的话可能会有用 popcount(A[]) popcount(A[] & B[]) 或者对于不同的微体系结构,其中ALU与负载吞吐量不同。

另一个变体, 使用单个AVX512VBMI(Cannonlake?)可以统计每个字节7个元素 vpermi2b ( _mm512_permutex2var_epi8 ) vpshufb

使用AVX512VBMI计数压缩-8(但不包括 AVX512VPOPCNTDQ vpermi2b 要计算低7,请按住shift键并屏蔽顶部位,然后将其相加。(单个位的popcount=该位)。

uint8_t 元素更容易有效地洗牌(因为有字节洗牌,如 vpshufb ),因此,如果您必须在飞行中进行转置,则可能值得考虑。或者在转换时只在飞行中打包到位?

这对于缓存占用空间/内存带宽来说也是一件大事 .8大小差异的因素可能意味着整行或整列只占用一级缓存的一部分,而不是溢出一级缓存。因此,它可以使缓存阻塞更容易/不那么重要。

1MiB on Skylake-AVX512 ),但将适用于许多核心Xeon CPU上的L3。但是,L3由所有核心(包括云环境中的其他虚拟机)竞争性共享,并且比L2慢得多。(像您将在HPC/cloud系统中运行的许多核心Xeon的单核内存带宽低于四核桌面,因为三级缓存的延迟更高,并发性没有增加(请参阅 "latency-bound platforms" section of this answer 。在Xeon上驱动相同数量的内存带宽需要更多的内核,即使总吞吐量更高。但是如果你能让每个核心都在其私有L2上工作,你会收获很多。)

:您已经安排了循环,因此需要将一次布尔运算减少到非零计数。这是一件好事。

使用8位整数0/1元素,最多可以执行255个操作 vpaddb vpsadbw against an all-zero vector to horizontally add the bytes in a vector into 64-bit integers .然后将蓄能器与 vpaddq , then horizontally sum it

VPSHUFB -基于位切片的popcount。(参见 http://wm.ite.pl/articles/sse-popcount.html 例如如果必须手动对其进行矢量化,则需要使用内部函数而不是asm来编写。)

您可以考虑将数据打包为每字节4位的低半字节。 vpshufb 可以计算AND结果的每个字节中的位,而不需要任何移位/掩蔽。在内环中,有两个负载, , vpaddb 。通过适当的展开,这应该与每个时钟2x 32B的L1D负载带宽保持一致,并使所有三个矢量执行端口(在Haswell或Skylake上)饱和。每隔128个或255个向量或其他什么来累积累加器的字节 vpsadbw 因此,最内部的循环应该以每字节4个元素*每向量32B=每时钟周期128个元素的速度运行, 如果您可以安排它读取L1D缓存中的热数据。预计Haswell/Skylake上的二级缓存的带宽将达到一半左右,而三级缓存的带宽则要差得多。

具有如果元素为0或1,则可以使用一些整数乘加指令。它们的设计有点奇怪,用于不同的用例,而不是FP FMA。它们将乘法结果的水平对相加,生成更宽的元素。 VPMADDUBSW vpsadbw vpsadbw ,这对你毫无益处 vpand 。只有当您想使用 vpaddw vpmaddubsw doesn't seem useful here, because

使用SSE/AVX可以有效地将0/1整数转换为位图 :对于32位整数元素, vpslld ymm0, 31 将相关位左移到每个元素的顶部,然后 vmovmskps eax, ymm0 获取每个32位元素的高字节的8位掩码。对于8位整数元素, vpslld ymm0, 7 / vpmovmskb eax, ymm0 执行相同的操作,但每个字节都会生成一个32位整数位图结果。(只有每个字节的符号位很重要,因此没有只有8位粒度的移位指令是很好的。您不需要对进入下一个元素的位做任何事情。)

这不是一种很好的立即使用向量的方法,因为最终结果会出现在整数寄存器中。这不是一种很好的即时生成和使用的格式,但它是最紧凑的,因此如果您可以长期保持这种格式的矩阵,那么它是有意义的。(如果加载时内存带宽有限。)

:单向是2x vpackssdw vpacksswb 。因为它们在128b车道内运行,您的元素最终将被重新排序。但只要每一行/列的顺序相同,就可以了。只有当您想要获取一个不是以32个元素的倍数开始的行/列块时,这才是一个问题。这里的另一个选项是左移(8、16和24),和或向量一起。事实上

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.
}

转换为每字节4位半压缩 :我们可以使用与上述相同的想法。从中获取4个向量 load_interleave4x32 uint8\u t set1_epi8(1)

如果以这种格式存储整个矩阵,则可以将其用作转置的一部分,或者可以使用这种格式存储缓存阻塞转置的临时副本。matmul会多次接触每一行/每一列,因此第一次制作紧凑格式可能值得做额外的工作,因为这样可以在后续过程中为每个向量做4倍的工作。

带AVX512BW(Skylake-AVX512)

popcnt

@Harold points out an interesting identity 这让我们能够以额外的整数运算为代价,进行2/3倍多的向量popcounts。

   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 2* 如果我们使用单独的预移位 vpshufb LUT向量。另请参见 this implementation that uses 30x vpternlogd + 1 vector popcnt to handle 16 vectors of 512b 16*popcnt 循环内计数;其他一切都是连锁的)。

这对于计算完全压缩的每字节8位元素很可能是值得的,并且使该格式对AVX512更有吸引力,相比之下,为没有太多移位/掩蔽的popcounting优化的密度较低的格式。

VPERNLOGD 如果是字节粒度,则也可以用作转置的位混合指令 VPBLENDMB zmm{k1}, zmm, zmm 纹理不够精细。

对于某些CPU上的AVX2来说,这可能是值得的,也许可以避免每4或5个向量popcounts中就有1个,而不是3个向量popcounts中的1个?或者,如果它只是增加了执行端口的总压力,并且在任何特定端口上都没有瓶颈,那么它可能根本没有帮助。这对标量很有用指令(可能在没有AVX2的CPU上),因为这些指令在Intel CPU的单个端口上会出现瓶颈。

uint8\u t 将布尔元素转换为非交错位图的效率略高于AVX2(甚至不需要移位),相反的转换效率更高。测试掩码或将掩码与set1\u epi8(1)的向量进行比较都可以完成这项工作,从64字节的输入中生成64位掩码。或者从32位整数开始,一次生成16位掩码。您可以有效地将这些位与 kunpck

_mm512_test_epi8_mask ( vptestmb ) 有趣的是:将两个向量合并在一起,并生成真/假字节元素的掩码寄存器结果。但这并不是我们真正想要的:如果我们要打包我们的位,我们希望将其作为输入矩阵的预处理步骤,而不是在进行内积时即时进行。

位图->0/-1的向量很快: __m512i _mm512_movm_epi8 (__mmask64 k) ( vpmovm2b ) -1 而不是添加 1 ,但必须先屏蔽它,然后才能将一个字节中的多个位合并在一起。

,您没有512b vpshufb AVX512 popcnt extension 直接用于矢量popcnt,但尚未宣布其硬件。(AVX2 vpshufb ymm 但是,KNL速度非常慢:每12个周期一次,并且 psadbw ymm a bithack popcnt based on 32-bit integer elements, since that's just AND/shift/ADD .32位元素执行popcnt的步骤将少于64位,但仍然足够大,不会溢出合理的问题大小(因此可以将向量的水平和延迟到循环外)

考虑到存储格式的选择,对KNL来说,每字节打包多个位可能不是一个好主意,但单字节整数元素很好。 vpandd zmm 和 vpaddd zmm SWAR

转置位:

BMI2 _pdep_u64

另一个有用的选择是 vpmovmskb this blog post How would you transpose a binary matrix? .

在matmul中使用

有些选择取决于输入数据的格式,以及重复使用相同矩阵的频率。如果一个矩阵将被多次使用,那么提前将其压缩到每字节4或8位是有意义的。(或第一次使用时正在飞行中)。保留它的转置副本也可能有意义,尤其是当它始终是需要转置的乘法的一侧时。(如果您有时需要一种方式,有时需要另一种方式,那么动态重做可能会更好地占用三级缓存空间。但这些缓存空间足够大,您可能不会获得太多三级命中率,因此只保留一个转置副本就可以了。)

或者甚至可以在从输入格式转换时写出转置和非转置版本。

当缓存阻塞一个普通的FP matmul时,同样的原则也适用,所以请阅读这篇文章。

使用位集 & 因为整个列将把值放回内存中,然后在中再次循环 .count() 关于结果。我怀疑编译器是否会将其优化为使用 VPSHUFB -基于每个矢量的位切片popcnt VPAND http://wm.ite.pl/articles/sse-popcount.html

根据矩阵大小,至少内部循环可能会命中L1D缓存,但两次循环产生的额外加载/存储指令开销更大,并且还会干扰有价值数据的预取。

这并不容易。唯一不糟糕的是 clang++ -stdlib=libc++ vector<bool> std::count(v.begin(), v.end(), true); 到矢量化 vpshufb + + vpaddb 在展开的环内 vpsadbw + vpaddq 每次迭代一次,但这对于自动矢量化的代码非常好。

向量(<);布尔> 也是位图,但是 std::count(v.begin(),v.end(),true); 非常糟糕:它使用了一个完全幼稚的循环,每次测试1位。它甚至不能有效地做到这一点。相同的 clang++ 默认情况下 libstdc++ 而不是新的 libc++ .

boost::dynamic_bitset .count() 成员函数,但它没有利用 popcnt公司 指令或AVX2。它执行逐字节LUT查找。这比 std::count(vector<bool>) 没有libc++,但对于HPC来说,它甚至不够好。

这是测试代码 on the Godbolt compiler explorer ,具有gcc和clang asm输出。他们都用过 -march=haswell .

但不幸的是,似乎没有一种有效的方法来实现按位和2 std::vector<bool> This answer 展示了如何获得g++的底层实现 ,但该代码不会自动矢量化。做同样的事情并对其进行调整,使其自动矢量化可以 向量(<);布尔> ,因为向量的向量是一个糟糕的额外间接级别。如果问题的转置部分也对性能至关重要,那么使用标准容器访问有效的popcount并不能解决整个问题。

对于 std::bitset<1024*1024>.count() 锂离子电池++ popcnt公司 说明,其中(根据 this

On vector<bool> â Howard Hinnant ,了解有关C++标准库的一些评论,以及为什么位数组是一种有用的数据结构,但是 向量(<);布尔> bool -每字节 bool[] 数组,与朴素的 (就像没有libc++的gcc和clang一样)。