特征:如何防止大型对象的额外副本;在RHS上未实现完整矩阵的情况下分配结果
Eigen: how to prevent extra copies of a large object; assign to result without realizing full matrix on RHS
如果这是一些我没能理解的基本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;
}
您应该重写表达式。从数学上讲,从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版本(还有许多其他改进)。
我无法编译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')
- 为什么 MSVC _count_of实现将 0 添加到 sizeof 的结果中?
- OpenCV 混合模式实现:为什么看似等效的操作会产生不同的结果?
- libc++ 对 std::map/set::equal_range 的实现给出了意想不到的结果
- 微小加密算法实现会产生意想不到的结果
- 为什么在 C++ 中实现高斯勒让德算法没有产生结果
- 反向 Cuthill-McKee 算法:在 Matlab 和提升实现中具有无与伦比的结果
- 在我的C++链表实现中取消引用节点指针,给出意想不到的结果
- 在不同的平台/编译器上实现相同的浮点计算结果
- 为什么定义了sizeof实现的结果
- 特质实现给出了Clang和G 的不同结果,这是正确的
- C void*任何类型实现都会返回奇怪的结果
- Negamax C++实现给出了错误的结果
- 着色器中实现双精度cos()的结果是NaN,但在CPU上运行良好.出了什么问题
- C++,Lapack-Cholesky分解实现结果不准确
- 如何在c++中实现此结果?指向数组的数组
- 实现多个返回语句 c++ 时出现奇怪的结果
- 绕过虚拟模板功能以实现期望的结果
- C++成员函数结果缓存实现
- C++中合并排序实现的特殊行为,排序后的结果不应该在适当的位置
- HashTable实现的错误结果