在 C++ 中处理大型矩阵答案

【问题标题】：Handling large matrices in C++在 C++ 中处理大型矩阵
【发布时间】：2015-06-06 18:47:25
【问题描述】：

我在 C++ 中使用 doubles 的大型矩阵。我需要从这些矩阵中获取行或列并将它们传递给函数。我能做到这一点的最快方法是什么？

一种方法是编写一个函数，将所需行或列的副本作为std::vector 返回。
另一种方法是将整个事物作为参考传递并修改函数以能够读取所需的值。

还有其他选择吗？你推荐哪一个？

顺便说一句，你建议我如何将数据存储在矩阵类中？我现在正在使用std::vector< std::vector< double > >。

编辑

我必须提到矩阵可能有两个以上的维度。所以在这里使用 boost 或 arma::mat 是不可能的。虽然，我在图书馆的其他部分使用犰狳。

【问题讨论】：

vector 是一个糟糕的选择 - 只需使用单个向量并提供运算符E operator () (row, column)
如果您的矩阵是vector<vector<double>>，那么您可以传递对行的引用（即vector<double>&）
我还没有尝试过，但是 boost 1.58 有一个线性代数库，它有一个密集和稀疏的矩阵实现
这真的取决于函数在做什么以及矩阵有多大。例如，如果函数多次访问列的元素，那么制作副本可能是您的最佳选择，因为它会更好地使用缓存。在性能方面，您希望衡量您的实际用例。
犰狳有cube class，可以像“3D矩阵”一样对待

标签： c++ c++11 matrix

【解决方案1】：

如果大于 2 的可变维数是关键要求，请查看boost's multidimensional array library。它具有高效（免费复制）“views”，您可以使用它来引用整个矩阵的低维“切片”。

这类事情“最快”的细节在很大程度上取决于您到底在做什么，以及访问模式/工作集“足迹”如何适合您的硬件的不同级别的缓存和内存延迟；在实践中，值得复制到更紧凑的表示以获得更多的缓存一致访问，而不是进行稀疏跨步访问，这只会浪费大量缓存行。替代方案是Morton-order 访问方案，它至少可以在所有轴上摊销“坏轴”效应。不过，只有您自己对自己的代码和硬件用例进行基准测试才能真正回答这个问题。

（请注意，我不会将 Boost.MultiArray 用于 2 维数组 - 对于线性代数/图像处理应用程序有更快、更好的选择 - 但对于 3+ 来说，值得考虑。）

【讨论】：

【解决方案2】：

我会使用像http://arma.sourceforge.net/ 这样的库，因为您不仅可以获得存储矩阵的方法。你也有可以对其进行操作的函数。

【讨论】：

Armadillo（上面提到的库）是一个不错的选择——尤其是当您需要与 MATLAB 交互时。 Eigen3 也很强大，而且只有 header，这很方便。

【解决方案3】：

高效（多）线性代数是一门非常深奥的学科；没有简单的万能答案。主要挑战是数据局部性：您的计算机的内存硬件已针对访问连续的内存区域进行了优化，并且可能一次无法在缓存行之外的任何东西上运行（即使可以，效率也会下降）。

高速缓存行的大小各不相同，但可以考虑 64 或 128 字节。

因此，将数据布置在矩阵中以便可以在多个方向上高效地访问是一项不小的挑战；对于更高等级的张量更是如此。

此外，最佳选择可能在很大程度上取决于您对矩阵所做的确切操作。

你的问题真的不是像这样的问答形式能得到满意回答的。

但至少要让您开始研究，这里有两个关键词可能值得研究：

块矩阵
快速转置算法

您最好使用库而不是尝试自己创建库；例如闪电战++。（免责声明：我没有使用过blitz++）

【讨论】：

【解决方案4】：

vector<vector<...>> 分配速度慢，释放速度慢，访问速度慢，因为它有多个取消引用（不适合缓存）。

仅当您的（行或列）大小不同时（jagged arrays），我才会推荐它。

对于“正常”矩阵，您可以选择以下内容：

template <class T, size_t nDim> struct tensor {
    size_t dims[nDim];
    vector<T> vect;
};

并重载operator(size_t i, size_t j, etc.) 以访问元素。

operator() 必须进行索引计算（您必须在row-major 或列优先顺序之间进行选择）。对于nDim > 2，它变得有些复杂，它可以从缓存一些索引计算中受益。

要返回一行或一列，您可以定义子类型。

template <class T, size_t nDim> struct row /*or column*/ {
    tensor<T, nDim> & tensor;
    size_t iStart;
    size_t stride;
}

然后定义一个operator(size_t i)，它将返回tensor.vect[iStart + i*stride]

stride 值将取决于它是行还是列，以及您的（行优先或列优先）排序选择。

stride 将是 1 对于您的子类型之一。请注意，对于这种子类型，迭代可能会快得多，因为它对缓存友好。不幸的是，对于其他子类型，它可能会相当慢，而且您无能为力。

请参阅其他 SO 问题，了解为什么先迭代行然后列可能会比先迭代列然后再迭代具有巨大的性能差异。

【讨论】：

【解决方案5】：

我建议您通过引用传递它，因为复制可能是一个缓慢的过程，具体取决于大小。如果您想要扩展和收缩容器的能力，则 std::vector 很好。

【讨论】：