特征:如何防止大型对象的额外副本;在RHS上未实现完整矩阵的情况下分配结果

Question

如果这是一些我没能理解的基本c++，我提前道歉。

在展示我的代码之前，让我解释一下我想要完成什么。我有一个稀疏矩阵U和一个向量r，我想计算(U-r)(U-r)'其中减去U的每一列

然而，我不能一次完成所有这些，因为U-r是密集的，并且会增加内存使用(~ 700万列vs ~ 20,000行)。

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

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

因为我想尽可能地并行化这个操作(这是相当昂贵的)，所以减少内存使用是至关重要的。

那么，第一步是让我从3个拷贝到2个拷贝。

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

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

最后是代码。我在R里做并行化，这段代码只代表了并行化的一个过程。我正在传递一个列索引的连续向量列表来计算。

// [[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具有相同列数的行向量)。展开给你:

(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())，并存储为致密向量。剩下的操作是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);

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

Answer 2

我无法编译chtz的代码。我本想给他们的答案加分，但用户Michael Albers认为编辑回复以包含正确的代码是不可接受的。所以我必须用正确的答案创建一个新的帖子。

在转换成密集矩阵之前，我必须为U的外积创建一个中间稀疏矩阵。这似乎不太理想，我看到其他人有这个问题，但没有办法解决它。在任何情况下，这个结果将编译为:

// [[Rcpp::export]]
MatrixXd gen_Sigma_cpp_sp(const MSpMat &U, const MVec &r) {
    VectorXd UcolSum = U * VectorXd::Ones(U.cols());
    MatrixXd S = MatrixXd(SparseMatrix<double>(U * U.transpose())).
                    selfadjointView<Lower>().rankUpdate(UcolSum, r, -1.0).
                                             rankUpdate(r, U.cols());
    return S;
}

对于任何在R中使用这个的人，我必须在强制输入'dpoMatrix'之前将其包装在forceSymmetric中，这是一个普通的交叉prod(U - R)会给出的，并且对计算最有帮助:

SigmaS0 = as(forceSymmetric(gen_Sigma_cpp_sp(U, r), 'L'), 'dpoMatrix')