c++ - Eigen:如何防止大对象的额外拷贝；分配给结果而不实现 RHS 上的完整矩阵

Question

如果其中一些是我未能理解的基本 C++，我提前道歉。

在展示我的代码之前，让我解释一下我想要完成的事情。我有一个稀疏矩阵 U 和一个 vector r，我想计算 (U-r)(U-r)'，其中减法是针对 U 的每一列。

但是，我无法一次完成所有操作，因为 U-r 非常密集并且内存使用量激增(大约 700 万列与大约 20,000 行)。

利用外积 XX' 可以一次计算一列的事实，XX' == sum(XcXc')，其中 sum 是矩阵加法，我的策略是取几列，做减法和外积并累加结果。一次仅使用几列可将内存使用量降低到非常合理的数量(几百 MB)。

从表面上看，这需要 2 份 20,000 x 20,000 矩阵(每个 3.5 GB)，一份用于累积结果，一份用于临时右侧。但是，由于我不明白的原因，根据观察到的内存使用情况，我有 3 个拷贝。

因为我想尽可能多地并行化此操作(这非常昂贵)，所以减少内存使用量至关重要。

所以，第 1 步是让我从 3 份到 2 份。

如果可能的话，第 2 步是要认识到没有理由永远不需要在 RHS 上实现结果。也就是说，没有理由不继续将计算结果添加到累加矩阵的每个元素，而不是在 RHS 上创建一个临时矩阵，然后执行累加器矩阵的加法。

第 3 步是利用生成对称矩阵这一事实来减少计算时间。我认为这是通过 .selfadjointView(Lower) 完成的，但我无法准确解析如何在一致的基础上继续这样做。

最后是代码。我在 R 中进行并行化，这段代码只代表一个并行化过程。我正在传递要计算的列索引的连续 vector 列表。

// [[Rcpp::depends(RcppEigen)]] 
#include <iostream>
#include "Rcpp.h"
#include "RcppEigen.h"
#include "Eigen/Dense"
#include "Eigen/Sparse"

using Eigen::MatrixXd;

typedef Eigen::MappedSparseMatrix<double> MSpMat;
typedef Eigen::Map<Eigen::VectorXd> MVec;
typedef Eigen::Map<MatrixXd> MMat;



/*
 * tcrossprod_cpp just compute X * X' where X is a matrix, * is matrix
 * multiplication and ' is transpose, but in an efficient manner,
 * although it appears that R's tcrossprod is actually faster. Pulled it from
 * the RcppEigen book.
 */


MatrixXd tcrossprod_cpp(const MatrixXd &U) {
    const long m(U.rows());
    MatrixXd UUt(MatrixXd(m, m).setZero().
            selfadjointView<Eigen::Lower>().rankUpdate(U));
    return UUt;
}

// [[Rcpp::export]]
MatrixXd gen_Sigma_cpp_block_sp(const Rcpp::List &col_list, const MSpMat &U,
                                const MVec &r, int index1 = 1) {
    long nrow = U.rows();
    MatrixXd out = MatrixXd::Constant(nrow, nrow, 0.0);
    long ncol;
    Rcpp::IntegerVector y;
    for (long i = 0; i < col_list.size(); i++) {
        if (i % 10 == 0) {
            Rcpp::checkUserInterrupt();
        }
        y = col_list[i];
        ncol = y[y.size() - 1] - y[0] + 1;
        out.noalias() += tcrossprod_cpp((MatrixXd (U.block(0, y[0] - index1,
                                         nrow, ncol))).colwise() - r);
    }
    return out;
}

Answer 1

你应该重写你的表达。从数学上讲，从 U 的每一列中减去 r 与 U - r*ones 相同(其中 ones 是列数与 U 相同的行 vector )。扩展为您提供:

(U-r*ones)*(U-r*ones)^T = U*U^T - (U*ones^T)*r^T - r*(ones*U^T) + r*(ones*ones^T)*r^T

ones*ones^T等于U.cols()，U*ones^T可以计算为 U*VectorXd::Ones(U.cols()) 并存储到一个密集 vector 中。剩余的操作是 U*U.transpose() 的一个稀疏乘积(您可以直接将其存储到密集矩阵中，因为您的最终结果将是密集的，然后是两个等级更新:

VectorXd Usum = U * VectorXd::Ones(U.cols()); // sum of columns of U
MatrixXd result = U*U.transpose();
result.selfadjointView<Lower>().rankUpdate(Usum, r, -1.0);
result.selfadjointView<Lower>().rankUpdate(r,U.cols());

回答关于额外临时工的问题:在 tcrossprod_cpp 中，您创建了一个临时的 MatrixXd(m,m) 并将结果存储到 MatrixXd UUt 中。你其实完全可以避免这种方法，直接写

out.selfadjointView<Lower>().rankUpdate(MatrixXd(U.block(0, y[0] - index1,
                                     nrow, ncol))).colwise() - r);

编辑:在 Eigen 3.3(我正在测试 3.3rc1)之前，将稀疏乘积直接分配给密集矩阵显然是不可能的。如果可能的话，我建议切换到 3.3 版(还有许多其他改进)。