稀疏矩阵-密集向量与编译时已知矩阵的乘法
Sparse matrix-dense vector multiplication with matrix known at compile time
我有一个稀疏矩阵，其中只有零和一个作为条目（例如，形状为 32k x 64k 和 0.01% 的非零条目，并且没有模式可以利用非零条目的位置).该矩阵在编译时已知。我想用包含 50% 1 和 0 的非稀疏向量（编译时未知）执行矩阵向量乘法（模 2）。我希望这是高效的，特别是，我正在尝试利用矩阵在编译时已知的事实。
以高效格式存储矩阵（仅保存“ones”的索引）总是需要几兆字节的内存，直接将矩阵嵌入到可执行文件中对我来说似乎是个好主意。我的第一个想法是自动生成 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++ 代码相同的二进制文件。 您认为这值得一试吗（我想我还是会在星期一试一试）？
另一个想法是尝试将输入（也可能是输出？）向量存储为打包位并执行类似的计算。我希望人们无法绕过大量的位移位或与操作，这最终会变得更慢、更糟糕。
对于如何做到这一点，您还有其他想法吗？
I'm not sure why though - is the sheer codesize slowing it down that much?
问题是可执行文件很大，OS 会从您的 存储设备 中获取大量页面。这个过程非常缓慢。处理器通常会停止等待数据加载。甚至代码已经加载到 RAM 中（OS 缓存），这将是低效的，因为 RAM 的速度（延迟 + 吞吐量）非常糟糕。这里的主要问题是所有指令只执行一次。如果你重用这个函数，那么代码需要从缓存中重新加载，如果它太大而无法放入缓存，它将从慢速 RAM 中加载。因此，与其实际执行相比，加载代码的开销非常高。要克服这个问题，您需要使用相当 小的代码，循环迭代相当少量的数据 。
Does this kind of code rule out compiler optimization?
这取决于编译器，但大多数主流编译器（例如 GCC 或 Clang）都会以相同的方式优化代码（因此编译时间较慢）。
Do you think this is worth trying (I guess I will try anyway on monday)?
是的，这个解决方案显然更好，尤其是当索引以紧凑的方式存储时。在您的情况下，您可以使用 uint16_t 类型存储它们。所有索引都可以放在一个大缓冲区中。每行索引的 starting/ending 位置可以在引用第一个（或使用指针）的另一个缓冲区中指定。该缓冲区可以在应用程序开始时从专用文件加载到内存中一次，以减小生成的程序的大小（并避免在关键循环中从存储设备中获取数据）。具有非零值的概率为 0.01%，生成的数据结构将占用不到 500 KiB 的 RAM。在普通的主流桌面处理器上，它可以放入 L3 缓存（相当快），我认为假设 multiply 的代码经过 仔细优化，您的计算时间不应超过 1ms .
Another idea would be to try storing the input (and maybe also output?) vector as packed bits and perform the calculation like that.
只有当您的矩阵不是太稀疏时，位打包才有用。对于一个填充了 50% 非零值的矩阵，位打包方法非常棒。对于 0.01% 的非零值，位打包方法显然很糟糕，因为它需要太多 space.
I would expect one can't get around a lot of bit-shifting or and-operations and this would end up being slower and worse overall.
如前所述，从存储设备或RAM加载数据非常慢。在任何现代主流处理器上进行一些位移都非常快（并且比加载数据快得多）。
以下是计算机可以执行的各种操作的大致时间：
我实现了第二种方法（constexpr 数组以压缩列存储格式存储矩阵），它好多了。用 -O3 编译需要（对于包含 35'000 个的 64'000 x 22'000 二进制矩阵）<1min 并在我的笔记本电脑上在 <300 微秒内执行一次乘法（Julia 大约需要 350 微秒）计算）。可执行文件的总大小约为 1 MB。
可能还可以做得更好。如果有人有想法，请告诉我！
下面是一个代码示例（显示一个 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,0x5,0x5,0x5,0x5
};
constexpr std::array<std::uint16_t, num_nz> row_idx = {
0x0,0x1,0x2,0x3,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, 0, 1, 0, 1, 1, 0, 1, 0, 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;
}
我有一个稀疏矩阵，其中只有零和一个作为条目（例如，形状为 32k x 64k 和 0.01% 的非零条目，并且没有模式可以利用非零条目的位置).该矩阵在编译时已知。我想用包含 50% 1 和 0 的非稀疏向量（编译时未知）执行矩阵向量乘法（模 2）。我希望这是高效的，特别是，我正在尝试利用矩阵在编译时已知的事实。
以高效格式存储矩阵（仅保存“ones”的索引）总是需要几兆字节的内存，直接将矩阵嵌入到可执行文件中对我来说似乎是个好主意。我的第一个想法是自动生成 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++ 代码相同的二进制文件。 您认为这值得一试吗（我想我还是会在星期一试一试）？
另一个想法是尝试将输入（也可能是输出？）向量存储为打包位并执行类似的计算。我希望人们无法绕过大量的位移位或与操作，这最终会变得更慢、更糟糕。
对于如何做到这一点，您还有其他想法吗？
I'm not sure why though - is the sheer codesize slowing it down that much?
问题是可执行文件很大，OS 会从您的 存储设备 中获取大量页面。这个过程非常缓慢。处理器通常会停止等待数据加载。甚至代码已经加载到 RAM 中（OS 缓存），这将是低效的，因为 RAM 的速度（延迟 + 吞吐量）非常糟糕。这里的主要问题是所有指令只执行一次。如果你重用这个函数，那么代码需要从缓存中重新加载，如果它太大而无法放入缓存，它将从慢速 RAM 中加载。因此，与其实际执行相比，加载代码的开销非常高。要克服这个问题，您需要使用相当 小的代码，循环迭代相当少量的数据 。
Does this kind of code rule out compiler optimization?
这取决于编译器，但大多数主流编译器（例如 GCC 或 Clang）都会以相同的方式优化代码（因此编译时间较慢）。
Do you think this is worth trying (I guess I will try anyway on monday)?
是的，这个解决方案显然更好，尤其是当索引以紧凑的方式存储时。在您的情况下，您可以使用 uint16_t 类型存储它们。所有索引都可以放在一个大缓冲区中。每行索引的 starting/ending 位置可以在引用第一个（或使用指针）的另一个缓冲区中指定。该缓冲区可以在应用程序开始时从专用文件加载到内存中一次，以减小生成的程序的大小（并避免在关键循环中从存储设备中获取数据）。具有非零值的概率为 0.01%，生成的数据结构将占用不到 500 KiB 的 RAM。在普通的主流桌面处理器上，它可以放入 L3 缓存（相当快），我认为假设 multiply 的代码经过 仔细优化，您的计算时间不应超过 1ms .
Another idea would be to try storing the input (and maybe also output?) vector as packed bits and perform the calculation like that.
只有当您的矩阵不是太稀疏时，位打包才有用。对于一个填充了 50% 非零值的矩阵，位打包方法非常棒。对于 0.01% 的非零值，位打包方法显然很糟糕，因为它需要太多 space.
I would expect one can't get around a lot of bit-shifting or and-operations and this would end up being slower and worse overall.
如前所述，从存储设备或RAM加载数据非常慢。在任何现代主流处理器上进行一些位移都非常快（并且比加载数据快得多）。
以下是计算机可以执行的各种操作的大致时间：
我实现了第二种方法（constexpr 数组以压缩列存储格式存储矩阵），它好多了。用 -O3 编译需要（对于包含 35'000 个的 64'000 x 22'000 二进制矩阵）<1min 并在我的笔记本电脑上在 <300 微秒内执行一次乘法（Julia 大约需要 350 微秒）计算）。可执行文件的总大小约为 1 MB。
可能还可以做得更好。如果有人有想法，请告诉我！
下面是一个代码示例（显示一个 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,0x5,0x5,0x5,0x5
};
constexpr std::array<std::uint16_t, num_nz> row_idx = {
0x0,0x1,0x2,0x3,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, 0, 1, 0, 1, 1, 0, 1, 0, 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;
}