快速整数矩阵乘法与 bit-twiddling hacks

Question

我问是否有可能 显着改善 整数矩阵乘法与 bitwise operations。矩阵很小，元素是小的非负整数（小意味着最多20）。

为了让我们集中注意力，让我们非常具体，假设我有两个 3x3 矩阵，整数项为 0<=x<15。

以下天真的 C++ 实现执行一百万次执行大约 1 秒，用 linux time.

测量

#include <random>

int main() {
//Random number generator
std::random_device rd;
std::mt19937 eng(rd());
std::uniform_int_distribution<> distr(0, 15);

int A[3][3];
int B[3][3];
int C[3][3];
for (int trials = 0; trials <= 1000000; trials++) {
    //Set up A[] and B[]
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            A[i][j] = distr(eng);
            B[i][j] = distr(eng);
            C[i][j] = 0;
        }
    }
    //Compute C[]=A[]*B[]
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            for (int k = 0; k < 3; ++k) {
                C[i][j] = C[i][j] + A[i][k] * B[k][j];
            }
        }
    }
}
return 0;
}

备注：

矩阵不一定是稀疏的。
Strassen-like 评论在这里没有帮助。
我们尽量不要使用环境观察，在这个 具体问题 中，矩阵 A[] 和 B[] 可以编码为单个 64 位整数。想想稍微大一点的矩阵会发生什么。
计算是单线程的。

相关：Binary matrix multiplication bit twiddling hack and What is the optimal algorithm for the game 2048?

Answer 1

如果您对大量矩阵执行此计算，您可能会发现减小数据大小可以显着提高性能：

#include <cstdint>
#include <cstdlib>

using T = std::uint_fast8_t;

void mpy(T A[3][3], T B[3][3], T C[3][3])
{
for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            for (int k = 0; k < 3; ++k) {
                C[i][j] = C[i][j] + A[i][k] * B[k][j];
            }
        }
    }
}

奔腾可以在一条指令中移动和符号扩展一个8位值。这意味着每个缓存行获得的矩阵数量是原来的 4 倍。

更新：好奇心被激起，我写了一个测试：

#include <random>
#include <utility>
#include <algorithm>
#include <chrono>
#include <iostream>
#include <typeinfo>

template<class T>
struct matrix
{
    static constexpr std::size_t rows = 3;
    static constexpr std::size_t cols = 3;
    static constexpr std::size_t size() { return rows * cols; }

    template<class Engine, class U>
    matrix(Engine& engine, std::uniform_int_distribution<U>& dist)
    : matrix(std::make_index_sequence<size()>(), engine, dist)
    {}

    template<class U>
    matrix(std::initializer_list<U> li)
    : matrix(std::make_index_sequence<size()>(), li)
    {

    }

    matrix()
    : _data { 0 }
    {}

    const T* operator[](std::size_t i) const {
        return std::addressof(_data[i * cols]);
    }

    T* operator[](std::size_t i) {
        return std::addressof(_data[i * cols]);
    }

private:

    template<std::size_t...Is, class U, class Engine>
    matrix(std::index_sequence<Is...>, Engine& eng, std::uniform_int_distribution<U>& dist)
    : _data { (void(Is), dist(eng))... }
    {}

    template<std::size_t...Is, class U>
    matrix(std::index_sequence<Is...>, std::initializer_list<U> li)
    : _data { ((Is < li.size()) ? *(li.begin() + Is) : 0)... }
    {}


    T _data[rows * cols];
};

template<class T>
matrix<T> operator*(const matrix<T>& A, const matrix<T>& B)
{
    matrix<T> C;
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            for (int k = 0; k < 3; ++k) {
                C[i][j] = C[i][j] + A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

static constexpr std::size_t test_size = 1000000;
template<class T, class Engine>
void fill(std::vector<matrix<T>>& v, Engine& eng, std::uniform_int_distribution<T>& dist)
{
    v.clear();
    v.reserve(test_size);
    generate_n(std::back_inserter(v), test_size,
               [&] { return matrix<T>(eng, dist); });
}

template<class T>
void test(std::random_device& rd)
{
    std::mt19937 eng(rd());
    std::uniform_int_distribution<T> distr(0, 15);

    std::vector<matrix<T>> As, Bs, Cs;
    fill(As, eng, distr);
    fill(Bs, eng, distr);
    fill(Cs, eng, distr);

    auto start = std::chrono::high_resolution_clock::now();
    auto ia = As.cbegin();
    auto ib = Bs.cbegin();
    for (auto&m : Cs)
    {
        m = *ia++ * *ib++;
    }
    auto stop = std::chrono::high_resolution_clock::now();

    auto diff = stop - start;
    auto millis = std::chrono::duration_cast<std::chrono::microseconds>(diff).count();
    std::cout << "for type " << typeid(T).name() << " time is " << millis << "us" << std::endl;

}

int main() {
    //Random number generator
    std::random_device rd;
    test<std::uint64_t>(rd);
    test<std::uint32_t>(rd);
    test<std::uint16_t>(rd);
    test<std::uint8_t>(rd);
}

示例输出（最近的 macbook pro，64 位，使用 -O3 编译）

for type y time is 32787us
for type j time is 15323us
for type t time is 14347us
for type h time is 31550us

总结：

在此平台上，int32 和 int16 被证明彼此一样快。 int64 和 int8 同样慢（8 位结果让我吃惊）。

结论：

一如既往，向编译器表达意图，让优化器做它的事情。如果程序在生产中运行太慢，请进行测量并优化最严重的问题。

Answer 2

您链接的问题是关于矩阵的，其中每个元素都是一位。对于一位值 a 和 b，a * b 完全等同于 a & b.

对于添加 2 位元素，从头开始添加可能是合理的（并且比解包更快），使用 XOR（无进位加法），然后使用 AND、移位和屏蔽进位生成进位元素边界。

当添加进位产生另一个进位时，第 3 位将需要检测。与使用 SIMD 相比，我不认为模拟 3 位加法器或乘法器会是一个胜利。如果没有 SIMD（即在带有 uint64_t 的纯 C 中），它可能有意义。对于加法，您可以尝试使用普通加法，然后尝试撤消元素边界之间的进位，而不是通过 XOR/AND/shift 操作自己构建加法器。

打包与解包到字节的存储格式

如果您有很多这样的小矩阵，将它们以压缩形式（例如打包的 4 位元素）存储在内存中可以帮助减少缓存占用空间/内存带宽。 4 位元素很容易解压到每个元素都在向量的单独字节元素中。

否则，每字节一个矩阵元素存储。从那里，如果需要，您可以轻松地将它们解压缩为每个元素 16 位或 32 位，具体取决于目标 SIMD 指令集提供的元素大小。您可以将一些矩阵以解压缩格式保留在局部变量中以在乘法运算中重用，但将它们打包回每个元素 4 位以存储在数组中。

编译器在 x86 的标量 C 代码中用 uint8_t 搞砸了。请参阅关于@Richard 回答的评论：gcc 和 clang 都喜欢将 mul r8 用于 uint8_t，这迫使它们将数据移动到 eax（隐含的 input/output 用于一个 -操作数相乘），而不是 .

uint8_t 版本实际上比 uint16_t 版本运行得更慢，尽管它的缓存占用空间只有一半。

您可能会从某种 SIMD 中获得最佳结果。

英特尔 SSSE3 有一个 vector byte multiply, but only with adding of adjacent elements。使用它需要将你的矩阵解包成一个向量，在行之间有一些零或其他东西，所以你不会从一行中获取数据与另一行中的数据混合在一起。幸运的是，pshufb 可以将元素归零以及复制它们。

更有用的是 SSE2 PMADDWD，如果您将每个矩阵元素解压缩到单独的 16 位向量元素中。因此，给定一个向量中的一行，以及另一个向量中的转置列，pmaddwd (_mm_madd_epi16) 是一个水平 add 远离为您提供 C[i][j].

您可以将多个 pmaddwd 结果打包到一个向量中，这样您就可以一次存储 C[i][0..2]，而不是分别执行这些添加。

快速整数矩阵乘法与 bit-twiddling hacks

Fast integer matrix multiplication with bit-twiddling hacks

c++

algorithm

performance

matrix-multiplication

打包与解包到字节的存储格式

您可能会从某种 SIMD 中获得最佳结果。