如何在 C++ 中存储 100000*100000 矩阵？答案

【问题标题】：How do I store a 100000*100000 matrix in C++?如何在 C++ 中存储 100000*100000 矩阵？
【发布时间】：2016-02-08 09:53:51
【问题描述】：

我有两个向量 a[100000] 和 b[100000]。我想将a[i]*b[j] 存储在100000 x 100000 矩阵M 中。如何在 C++ 中做到这一点？

【问题讨论】：

如果你有足够的内存使用std::vector of std::vector，否则你可以使用memory mapped file。
不，我的意思是假设我正在制作一个矩阵 T 阶 n*n，其中每个元素 A(i,j)=A(i)*B(j)
还有很多库可以以节省空间的方式处理sparse matrixes。
@Joachim Pileborg 不，我必须在用户运行时存储它
所以？内存映射文件仍然是存储大量运行时数据的有效方式（我猜这就是您的意思）。没有什么说一旦您处理完数据就必须保留该文件。稀疏矩阵或（如果你有足够的内存）向量的向量，甚至单个向量都是其他有效的可能性。

标签： c++

【解决方案1】：

您可以使用std::vector<std::vector<your_type>> 来存储结果。

int rows = 100000, cols = 100000;

std::vector<std::vector<double>> result;
result.resize(rows);

for(int i=0; i<rows; i++) {
        result[i].resize(cols);
}    

for(int i=0; i<rows; i++) {
    for(int j=0; j<cols; j++){
        result[i][j] = a[i] * b[j];
    }
}

或者您可以使用线性代数库，例如 Eigen（您可能用它编写代码的时间会更少），这肯定会更高效。

【讨论】：

这是非常低效的。它具有糟糕的内存局部性和 100,000 次不必要的动态分配。首选大小为rows*cols 的单个平面向量，并在接口层覆盖二维索引方案。
至少每个子向量中缺少.reserve()，您可以（并且应该！）很容易地修复。
“效率低下”大约是每 MB 数据 40 个字节。 (0.004%)。作为交换，您不需要 40 GB 的连续分配。
或者直接创建它：std::vector<std::vector<double>> result { rows, std::vector<double>(cols) };.
@MSalters：这不仅仅是所需的内存量效率低下。在“初始化”时存在不必要的性能问题，在使用时也存在不必要的性能问题，我已经解决了这两个问题。对我来说，这不属于“过早优化；忽略它”的范畴，而是属于“实际上没有理由使用这种幼稚且天生缓慢的技术”的范畴。

【解决方案2】：

应该重新研究此答案的 NON-contiguity 部分。它可能是错的。

如果您想处理大量元素，例如 100000*100000。我不建议使用vector of vectors，因为内部vectors 元素彼此之间存在“非连续”属性。小的push_back 可能会造成很多混乱。

我会使用带有包装器的单个 vector。更多信息请参见Clean ways to write multiple 'for' loops。

【讨论】：

为什么你认为非连续性不好？它在这里可能很有用。
mmmm 老实说，我并没有深思熟虑。这是我以前经常听到的。无论如何，我能想到一些原因。例如，您不能在一批中应用诸如 std::max_elemnt 之类的东西..您必须分别针对每一行。嗯，我正在重新考虑。这与连续性无关，因为您可以为 std::list 执行此操作，它与迭代器有关.. mmmm 再次，我不知道！
如果您可以将一起使用的元素保持在同一页面上（4KB 以内），则连续性会有所帮助。但在连续情况下，相邻行中的元素相隔 400 kB。如果你想要一个 40 GB 矩阵的最大元素，你最好检查std::thread::hardware_concurrency 你应该为此创建多少线程。然后在你的线程上划分 100.000 行，并在每一行上调用 std::max_element。
非常感谢。我将阅读有关您提到的方法的更多信息

【解决方案3】：

    #include <vector>
    class C
    {
    public:
        C(const std::vector<double>& a_, const std::vector<double>& b_)
            :a(a_),b(b_){};
        double operator()(size_t i, size_t j) const { return a[i]*b[j]; }
    private:
         std::vector<double> a, b;
    };

问题究竟是什么？

最初的问题询问将C(i,j)=A(i)*B(j) 保存到矩阵的方法。

从 OOP 的角度来看，这样的 matrix 可以定义为一个对象，其方法接受两个输入（i 和 j），并返回一个结果（ret=A(i)*B(j)）。

这可以使用嵌套数组订阅 (c[i][j])、线性数组索引 (c[i*100000+j]) 或函数 (c.get(i, j)) 来实现。第三种方式也可以简化为函子（c.operator()(i, j) 或c(i, j)）。

然后呢？

如果您同意以上所有内容，则三个接口中的任何一个都可以达到目的，或者至少部分达到目的（就像我在评论中提到的那样，如果矩阵只需要提供随机 read 访问它的元素）。然后我们继续实施其中之一，第三个是我的选择。

为什么要那样做？

我的观察是，计算返回值并不昂贵，那么为什么不在实际访问产品时“懒惰地”计算产品呢？

这样存储效率很高（内存使用从n^2减少到2n）。

在getter函数中隐藏乘法并不会显着增加访问时间（两次内存访问和一次乘法，与理想情况下只有一次内存访问相比，但两种情况都是常数时间，这种实现对缓存更加友好用于减少数据的大小）。

因此，不是保存产品，而是只保存输入，而是在访问特定元素时计算产品。

缺少什么？

虽然可以操作这个“矩阵”（通过更改成员a 和b），但它不允许将任意元素更改为任意值。

实现数组切片的成员函数（如c(0:10:end, 4)）也不存在，但可行。

测试代码

int main() {
    C c({1,2,3,4},{10,20,30,40}); // a={1,2,3,4}; b={10,20,30,40}
    cout << "3*30 "<<c(2,2);      // c(2, 2) = a[2]*b[2] = 3*30 = 90
    return 0;
}

演示

http://ideone.com/bZR7AU

【讨论】：

好吧，我认为这个答案值得投反对票，因为该类并不是真正的矩阵（它缺乏随机写访问，就像你不能添加/乘以另一个兼容大小的任意矩阵一样并将结果存储回来）。但是如果你只从这个矩阵中读取，那么我认为这不是一个糟糕的解决方案。恕我直言，这个问题本身就是一个 x-y 问题。
没有真正回答问题，缺乏解释
@Pierre 解释已添加。
我不确定答案。很可能计算 A[i,j] = B[i] * C[j] 只是更大任务中的一步。

【解决方案4】：

如果系统中的 RAM 少于 80GB（对于 100000×100000 的双精度矩阵），使用内存中的 std::vector<double> 可能不可行。

以下是使用 mmap 文件执行此操作的方法。请查看内联 cmets：

#include <sys/mman.h>
#include <stddef.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <stdio.h>

#define ROWS 1000
#define COLS 1000
#define FILENAME "./matrix.doubles"

int main(void)
{
    double (*matrix)[ROWS][COLS]; // pointer to our matrix
    int fd; // file descriptor of backing file

    // open backing file
    fd = open(FILENAME,
              O_CREAT | O_RDWR, // create (if absent) and/or read and writable
              S_IRUSR | S_IWUSR); // (only) user may read and write
    if (fd < 0) {
        perror("Could not open file");
        return 1;
    }

    if ((lseek(fd, sizeof(*matrix), SEEK_SET) == (off_t) -1) ||
        ftruncate(fd, sizeof(*matrix)) ||
        (lseek(fd, 0, SEEK_SET) == (off_t) -1)) {
        perror("Could not set file size.");
        return 1;
    }

    matrix = mmap(NULL, // I don't care were the address starts
                  sizeof(*matrix), // size of matrix in bytes
                  PROT_READ | PROT_WRITE, // readable and writable
                  MAP_PRIVATE, // we access the data exclusively
                  fd, // file descriptor of backing file
                  0); // offset
    if (matrix == MAP_FAILED) {
        perror("Could not mmap file.");
        return 1;
    }

    // operate on matrix
    for (unsigned row = 0; row < ROWS; ++row) {
        for (unsigned col = 0; col < COLS; ++col) {
            (*matrix)[row][col] = row * col;
        }
    }

    // close backing file
    munmap(matrix, sizeof(*matrix));
    close(fd);

    return 0;
}

这是纯 C 代码。您可以通过使用例如std::array<double, ROWS, COLS>&而不是裸数组，但我认为这个想法应该可以理解。

【讨论】：

【解决方案5】：

如果您可以即时计算 a[i]*b[j]，那么您应该这样做，原因有两个：

从一个巨大的矩阵中获取结果可能并不比动态计算两个双精度值的乘积更快。
10000x10000 双矩阵需要 80 GB 的存储空间（内存中或磁盘上），并且可能需要一些额外的工作来访问预先计算的数据。

在下面的示例中，如果我动态计算 N=20000 的两个双精度值的乘积，我会看到 30 倍的速度提升（使用 clang 3.8 在发布模式下编译）。

template <typename T>
void test_lookup(std::vector<T> &data, std::vector<size_t> &index,
                 std::vector<T> &results) {
    const size_t LOOP = index.size() / 2;
    for (size_t idx = 0; idx < LOOP; ++idx) {
        auto row = index[2 * idx];
        auto col = index[2 * idx + 1];
        results[idx] = data[col * LOOP + row];
    }
}

template <typename T>
void test_mul(std::vector<T> &x, std::vector<T> &y, std::vector<T> &results) {
    for (size_t idx = 0; idx < x.size(); ++idx) {
        results[idx] = x[idx] * y[idx];
    }
}

【讨论】：