稀疏矩阵-密集向量乘法与编译时已知的矩阵

问题描述

我有一个稀疏矩阵，其中只有 0 和 1 作为条目（例如，形状为 32k x 64k 和 0.01% 的非零条目，并且在非零条目的位置方面没有可利用的模式）。该矩阵在编译时已知。我想使用包含 50% 1 和 0 的非稀疏向量（在编译时未知）执行矩阵向量乘法（模 2）。我希望这很有效，特别是，我试图利用矩阵在编译时已知的事实。

以有效的格式存储矩阵（仅保存“1”的索引）总是需要几兆字节的内存，直接将矩阵嵌入到可执行文件中对我来说似乎是个好主意。我的第一个想法是自动生成 C++ 代码，将所有结果向量条目分配给正确输入条目的总和。这看起来像这样：

constexpr std::size_t N = 64'000;
constexpr std::size_t M = 32'000;

template<typename Bit>
void multiply(const std::array<Bit,N> &in,std::array<Bit,M> &out) {
    out[0] = (in[11200] + in[21960] + in[29430] + in[36850] + in[44352] + in[49019] + in[52014] + in[54585] + in[57077] + in[59238] + in[60360] + in[61120] + in[61867] + in[62608] + in[63352] ) % 2;
    out[1] = (in[1] + in[11201] + in[21961] + in[29431] + in[36851] + in[44353] + in[49020] + in[52015] + in[54586] + in[57078] + in[59239] + in[60361] + in[61121] + in[61868] + in[62609] + in[63353] ) % 2;
    out[2] = (in[11202] + in[21962] + in[29432] + in[36852] + in[44354] + in[49021] + in[52016] + in[54587] + in[57079] + in[59240] + in[60362] + in[61122] + in[61869] + in[62610] + in[63354] ) % 2;
    out[3] = (in[56836] + in[11203] + in[21963] + in[29433] + in[36853] + in[44355] + in[49022] + in[52017] + in[54588] + in[57080] + in[59241] + in[60110] + in[61123] + in[61870] + in[62588] + in[63355] ) % 2;
    // LOTS more of this...
    out[31999] = (in[10208] + in[21245] + in[29208] + in[36797] + in[40359] + in[48193] + in[52009] + in[54545] + in[56941] + in[59093] + in[60255] + in[61025] + in[61779] + in[62309] + in[62616] + in[63858] ) % 2;
}

这确实有效（需要很长时间才能编译）。然而，它实际上似乎很慢（比 Julia 中相同的稀疏向量矩阵乘法慢 10 倍以上），并且比我认为必要的要多得多地增加可执行文件的大小。我对 std::array 和 std::vector 都进行了尝试，并且各个条目（表示为 Bit）为 bool、std::uint8_t 和 int，没什么进展也罢。我还尝试用 XOR 替换模和加法。 总而言之，这是一个糟糕的主意。我不知道为什么 - 纯粹的代码大小是否会减慢它的速度？这种代码是否排除了编译器优化？

我还没有尝试过任何替代方案。我的下一个想法是将索引存储为编译时常量数组（仍然给我巨大的 .cpp 文件）并循环遍历它们。最初，我希望这样做会导致编译器优化生成与我自动生成的 C++ 代码相同的二进制文件。 你认为这值得尝试吗（我想我周一还是会尝试的）？

另一个想法是尝试将输入（也可能是输出？）向量存储为压缩位并执行这样的计算。我希望人们无法解决大量的位移或与运算，这最终会导致整体速度变慢且更糟。

对于如何做到这一点，您有其他想法吗？

解决方法

我不知道为什么 - 纯粹的代码大小是否会减慢它的速度？

问题在于可执行文件很大，操作系统将从您的存储设备中获取大量页面。这个过程非常缓慢。处理器通常会停止等待数据加载。即使代码已经加载到 RAM 中（操作系统缓存），它也会效率低下，因为 RAM 的速度（延迟 + 吞吐量）非常糟糕。这里的主要问题是所有指令只执行一次。如果您重用该函数，则需要从缓存中重新加载代码，如果它太大而无法放入缓存中，则会从慢速 RAM 中加载。因此，与实际执行相比，加载代码的开销非常高。为了克服这个问题，您需要使用相当带有循环迭代的相当少量数据的小代码。

这种代码是否排除了编译器优化？

这取决于编译器，但大多数主流编译器（例如 GCC 或 Clang）都会以相同的方式优化代码（因此编译时间很慢）。

你认为这值得尝试吗（我想我周一还是会尝试的）？

是的，这个解决方案显然更好，特别是如果索引以紧凑的方式存储。在您的情况下，您可以使用 uint16_t 类型存储它们。所有索引都可以放在一个大缓冲区中。每行索引的开始/结束位置可以在另一个引用第一个（或使用指针）的缓冲区中指定。该缓冲区可以在应用程序开始时从专用文件加载到内存中，以减少生成的程序的大小（并避免在关键循环中从存储设备获取）。具有非零值的概率为 0.01%，生成的数据结构将占用不到 500 KiB 的 RAM。在普通的主流桌面处理器上，它可以放入 L3 缓存（相当快），我认为假设 multiply 的代码经过仔细优化，您的计算时间不应超过 1 毫秒 em>.

另一个想法是尝试将输入（也可能是输出？）向量存储为压缩位并执行这样的计算。

只有当您的矩阵不是太稀疏时，位打包才是好的。矩阵填充了 50% 的非零值，位打包方法很棒。对于 0.01% 的非零值，位打包方法显然很糟糕，因为它占用了太多空间。

我希望人们无法绕过大量的位移或与运算，这最终会导致整体速度变慢且更糟。

如前所述，从存储设备或 RAM 加载数据非常慢。在任何现代主流处理器上进行一些位移都非常快（并且比加载数据快得多）。

以下是计算机可以执行的各种操作的大致时间：

,

我实现了第二种方法（constexpr 数组以压缩列存储格式存储矩阵），它要好得多。 （对于包含 35,000 个的 64'000 x 22'000 二进制矩阵）-O3 编译并在我的笔记本电脑上在

可能还是可以做得更好。如果有人有想法，请告诉我！

下面是一个代码示例（显示一个 5x10 矩阵），说明了我所做的。

#include <iostream>
#include <array>

// Compressed sparse column storage for binary matrix
constexpr std::size_t M = 5;
constexpr std::size_t N = 10;
constexpr std::size_t num_nz = 5;
constexpr std::array<std::uint16_t,N + 1> colptr = {
0x0,0x1,0x2,0x3,0x4,0x5,0x5
};
constexpr std::array<std::uint16_t,num_nz> row_idx = {
0x0,0x4
};

template<typename Bit>
constexpr void encode(const std::array<Bit,N>& in,std::array<Bit,M>& out) {

    for (std::size_t col = 0; col < N; col++) {
        for (std::size_t j = colptr[col]; j < colptr[col + 1]; j++) {
            out[row_idx[j]] = (static_cast<bool>(out[row_idx[j]]) != static_cast<bool>(in[col]));
        }
    }
}

int main() {
    using Bit = bool;
    std::array<Bit,N> input{1,1,1};
    std::array<Bit,M> output{};
    
    for (auto i : input) std::cout << i;
    std::cout << std::endl;

    encode(input,output);

    for (auto i : output) std::cout << i;
}