Eigen：如何防止大对象的额外副本；分配给结果而不实现 RHS 上的完整矩阵

Question

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

在展示我的代码之前，让我解释一下我想要完成的事情。我有一个稀疏矩阵 U 和一个向量 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 中进行并行化，这段代码只代表一个并行化过程。我正在传递要计算的列索引的连续向量列表。

// [[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);

编辑： 在 Eigen 3.3 之前显然不可能直接将稀疏乘积分配给密集矩阵（我正在使用 3.3rc1 进行测试）。如果可以，我建议切换到 3.3 版本（还有许多其他改进）。

Answer 2

我无法编译 chtz 的代码。我本想将他们的答案归功于他们，但用户 Michael Albers 决定编辑响应以包含正确的代码是不可接受的。所以我必须用正确的答案创建一个新的 post。

在转换为密集矩阵之前，我必须为 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 使用它的任何人，我必须先将其包装在 forceSymmetric 中，然后才能强制键入 'dpoMatrix'，这是普通的 tcrossprod(U - r) 会给出的，并有助于大多数计算都是在线的：

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

Eigen：如何防止大对象的额外副本；分配给结果而不实现 RHS 上的完整矩阵

Eigen: how to prevent extra copies of a large object; assign to result without realizing full matrix on RHS

c++

rcpp

eigen