使用 Eigen 和 C++ 做海量矩阵乘积的 colsum答案

【问题标题】：Using Eigen and C++ to do a colsum of massive matrix product使用 Eigen 和 C++ 做海量矩阵乘积的 colsum
【发布时间】：2015-06-05 19:15:41
【问题描述】：

我正在尝试计算 colsum(N * P)，其中 N 是稀疏的 1M x 2500 矩阵，P 是密集的 2500 x 1.5M 矩阵。我将 Eigen C++ 库与英特尔的 MKL 库一起使用。问题是矩阵 N*P 实际上不能存在于内存中，它太大了（~10 TB）。我的问题是，Eigen 是否能够通过惰性评估和并行性的某种组合来处理这种计算？它在这里说 Eigen 不会不必要地制作临时矩阵：http://eigen.tuxfamily.org/dox-devel/TopicLazyEvaluation.html

但是，Eigen 是否知道将 N * P 分段计算为实际适合内存的块？ IE：它必须执行 colsum(N * P_1) ++ colsum(N * P_2) ++ .. ++ colsum(N * P_n) 之类的操作，其中 P 按列分为 n 个不同的子矩阵，“+ +" 是连接。

我正在使用 128 GB RAM。

【问题讨论】：

标签： c++ matrix eigen

【解决方案1】：

我试了一下，但最终得到了一个错误的 malloc（我只在 Win8 上运行 8GB）。我设置了我的main() 并使用了我编写的非内联colsum 函数。

int main(int argc, char *argv[])
{
    Eigen::MatrixXd dense = Eigen::MatrixXd::Random(1000, 100000);
    Eigen::SparseMatrix<double> sparse(100000, 1000);
    typedef Triplet<int> Trip;
    std::vector<Trip> trps(dense.rows());
    for(int i = 0; i < dense.rows(); i++)
    {
        trps[i] = Trip(20*i, i, 2);
    }
    sparse.setFromTriplets(trps.begin(), trps.end());

    VectorXd res = colsum(sparse, dense);
    std::cout << res;
    std::cin >> argc;
    return 0;
}

尝试很简单：

__declspec(noinline) VectorXd 
colsum(const Eigen::SparseMatrix<double> &sparse, const Eigen::MatrixXd &dense)
{
    return (sparse * dense).colwise().sum();
}

这有一个糟糕的malloc。看来您必须自己手动拆分它（除非其他人有更好的解决方案）。

编辑

我稍微改进了功能，但得到了同样糟糕的 malloc：

__declspec(noinline) VectorXd 
colsum(const Eigen::SparseMatrix<double> &sparse, const Eigen::MatrixXd &dense)
{
    return (sparse * dense).topRows(4).colwise().sum();
}

编辑 2

另一种选择是使稀疏矩阵密集并强制进行惰性求值。我认为它不适用于稀疏矩阵（哦，好吧）。

__declspec(noinline) VectorXd 
colsum(const Eigen::SparseMatrix<double> &sparse, const Eigen::MatrixXd &dense)
{
    Eigen::MatrixXd denseSparse(sparse);
    return denseSparse.lazyProduct(dense).colwise().sum();
}

这不会给我糟糕的 malloc，但会计算出很多毫无意义的 0*x_i 表达式。

【讨论】：

【解决方案2】：

回答您的问题：特别是在涉及产品时，Eigen 经常将表达式的一部分评估为临时变量。在某些情况下，这可以优化但尚未实施，在某些情况下，这基本上是实现它的最有效方式。

但是，在您的情况下，您可以简单地计算 N（1 x 2500 向量）的 colsum 并将其乘以 P。也许未来版本的 Eigen 将能够自己进行这种优化，但大多数情况下，在让计算机完成其余工作之前，自己进行针对特定问题的优化是个好主意.

顺便说一句：恐怕sparse.colwise() 还没有实现，所以你必须手动计算。如果你很懒，你可以计算 Eigen::RowVectorXd Nsum = Eigen::RowVectorXd::Ones(N.rows())*P;（我没有检查过，但是这实际上可能会优化到接近最优的代码，使用最新版本的 Eigen）。

【讨论】：