一种快速的位旋转运算整数矩阵乘法方法

Question

我是在问是否有可能用按位运算大大改进整数矩阵乘法。矩阵是小的，元素是小的非负整数(小均值最多为20)。

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

下面这个简单的C++实现执行了一百万次，用linux time来衡量，大约执行了1s。

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

备注：

1. 矩阵不一定是稀疏的。
2. 斯特拉森样注释在这里没有帮助。
3. 让我们不要使用间接观察，在这个特定的问题中，矩阵A[]和B[]可以编码为单个 64位整数。想想只会发生在更大的矩阵上的事情。
4. 计算是单线程的.

​

相关：二进制矩阵乘法位旋转黑客和游戏2048的最优算法是什么？

Answer 1

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

对于添加2位元素，使用XOR (无进位添加)从零开始添加可能是合理的(而且比解压更快)，然后生成带进位和移位，并屏蔽跨元素边界的进位。

当添加进位产生另一个进位时，第三位将需要检测。我不认为这将是一个胜利，甚至模拟一个3位加法器或乘法器，比使用SIMD。没有SIMD (即在纯C和uint64_t)，这可能是有意义的。对于add，您可以尝试使用普通的add，然后尝试撤消元素边界之间的进位，而不是用XOR/ and /shift操作自己构建加法器。

打包和解压缩到字节存储格式

如果您有很多这些小矩阵，将它们以压缩的形式存储在内存中(例如，打包的4位元素)可以帮助提高缓存占用/内存带宽。4位元素很容易解压缩，以便将每个元素放在向量的单独字节元素中。

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

编译器在的标量C代码中使用了。参见关于@Richard的答案的评论: gcc和clang都喜欢使用mul r8 for uint8_t，这迫使他们将数据移动到eax (一个操作数的隐式输入/输出)，而不是忽略目标寄存器的低8位之外的垃圾。。

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

你可能会从某种SIMD中得到最好的结果。

英特尔SSSE3有一个向量字节相乘，但只需添加相邻元素。使用它将需要将矩阵解压到一个向量中，在行或其他之间有一些零，这样就不会从一行中获取数据，而不会从另一行中获取数据。幸运的是，pshufb不仅可以复制元素，还可以实现零元素。

如果您在一个单独的16位向量元素中解压缩到每个矩阵元素，那么PMADDWD更有可能是有用的。因此，给定一个向量中的一行，另一个向量中的一个转置列，pmaddwd (_mm_madd_epi16)是一个水平add，与给出C[i][j]所需的点乘积结果相去甚远。

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

Answer 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位值.这意味着每个缓存行得到的matricies是原来的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，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位的结果让我感到惊讶)。

结论：

和以往一样，向编译器表示意图，让优化器完成它的任务。如果程序在生产过程中运行得太慢，那么就对最坏的人进行测量和优化。