用位快速整数矩阵乘法

2018-06-14 02:44:08

我在问是否有可能通过按位运算来改善相当大的整数矩阵乘法。矩阵很小，元素是小的非负整数（小的意思是最多20个）。

为了保持我们的专注，让我们非常具体，并说我有两个3x3矩阵，整数条目0 <= x <15。

以下天真的C ++实现执行了一百万次，执行时间大约为1s，用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;
}

笔记：

矩阵不一定是稀疏的。

斯特拉森式的评论在这里没有帮助。

让我们尝试不使用间接观察，在这一特定问题的矩阵A[]和B[]可以被编码为一个 64位的整数。想想只有更大的矩阵会发生什么。

计算是单线程的。

相关：二进制矩阵乘法位twiddling hack和什么是2048年游戏的最优算法？

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

为了添加2位元素，它可能似乎是合理的（并且比拆包更快），从基本上从头开始添加XOR（carryless-add），然后通过元素边界生成带有AND，shift和mask的进位。

当添加进位产生另一进位时，第三位需要检测。与使用SIMD相比，我认为即使是3位加法器或乘法器也不会赢。没有SIMD（即在纯C中使用uint64_t ）它可能是有道理的。对于添加，您可以尝试使用正常添加，然后尝试撤消元素边界之间的进位，而不是通过异或/移位操作自行构建加法器。

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

如果你有很多这种微型矩阵，以压缩形式（例如，打包的4bit元素）将它们存储在内存中可以帮助缓存占用空间/内存带宽。 4bit元素非常容易解包，从而将每个元素放在一个单独的矢量字节元素中。

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

编译器在uint8_t中使用x86的标量C代码 。请参阅@ Richard的回答：gcc和clang都喜欢使用mul r8作为uint8_t ，这迫使他们将数据移入eax （单操作数乘法的隐式输入/输出），而不是使用imul r32, r32并忽略在目标寄存器的低8位之外的垃圾。

uint8_t版本实际上比uint16_t版本运行速度慢，即使它具有缓存足迹的一半。

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

英特尔SSSE3具有向量字节乘法，但只能添加相邻元素。使用它需要将你的矩阵拆分成一个向量，在行之间有一些零，或者别的东西，所以你不会从一行中得到与另一行数据混合的数据。幸运的是， pshufb可以将元素归零并复制它们。

如果您在单独的16位矢量元素中解开每个矩阵元素， PMADDWD可能有用的是SSE2 PMADDWD 。因此，给定一个向量中的一行和另一个向量中的转置列， pmaddwd （ _mm_madd_epi16 ）是向您提供C[i][j]所需的点积结果的一个水平add 。

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

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

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

Pentium可以在一条指令中移动并签名扩展一个8位值。这意味着你每个缓存行获得4倍的matricies。

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

#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位结果令我感到惊讶）。

结论：

与往常一样，向编译器表达意图并让优化器完成它的工作。如果程序在生产中运行速度过慢，请进行测量并优化最差的违规者。

链接地址: http://www.djcxy.com/p/40185.html

上一篇: Fast integer matrix multiplication with bit

下一篇: sided items based on the similarity of consecutive items