RcppEigen中的MappedSparseMatrix [英] MappedSparseMatrix in RcppEigen
问题描述
我想使用RcppEigen软件包中实现的共轭梯度算法来求解大型稀疏矩阵.
I want to use conjugate gradient algorithm implemented in RcppEigen package for solving large sparse matrices.
由于我是Rcpp和C ++的新手,所以我只是从密集矩阵开始的.
Since I am new to Rcpp and C++, I just started with the dense matrices.
// [[Rcpp::depends(RcppEigen)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <Eigen/IterativeLinearSolvers>
using Eigen::SparseMatrix;
using Eigen::MappedSparseMatrix;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;
using Rcpp::as;
using Eigen::ConjugateGradient;
typedef Eigen::MappedSparseMatrix<double> MSpMat;
// [[Rcpp::export]]
VectorXd getEigenValues(SEXP As, SEXP bs) {
const Map<MatrixXd> A(as<Map<MatrixXd> > (As));
const Map<VectorXd> b(as<Map<VectorXd> > (bs));
ConjugateGradient<MatrixXd> cg;
cg.compute(A);
VectorXd x=cg.solve(b);
return x;
}
这按预期工作.因此,我想将其扩展为适合稀疏矩阵.
This works as expected. Therefore, I thought to extend this to suit to sparse matrices.
// [[Rcpp::depends(RcppEigen)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <Eigen/IterativeLinearSolvers>
using Eigen::SparseMatrix;
using Eigen::MappedSparseMatrix;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;
using Rcpp::as;
using Eigen::ConjugateGradient;
typedef Eigen::MappedSparseMatrix<double> MSpMat;
// [[Rcpp::export]]
VectorXd getEigenValues(SEXP As, SEXP bs) {
const MSpMat A = as<MSpMat>(As);
const Map<VectorXd> b(as<Map<VectorXd> > (bs));
ConjugateGradient<SparseMatrix<double> > cg;
cg.compute(A);
VectorXd x=cg.solve(b);
return x;
}
但是,这往往会产生非常奇怪的结果.代码本身也没有给出任何错误.希望有人可以帮助我纠正此错误.
However, this tends to give really strange results. The code in itself is also not giving any errors. Hope someone could help me to correct this error.
谢谢.
推荐答案
您需要在Eigen :: ConjugateGradient函数中使用Eigen :: MappedSparseMatrix类型.尝试以下代码:
You need to use your Eigen::MappedSparseMatrix type in the Eigen::ConjugateGradient function. Try the following code:
#include <RcppEigen.h>
typedef Eigen::MappedSparseMatrix< double > mappedSparseMatrix ;
typedef Eigen::Map< Eigen::VectorXd > mappedVector ;
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Eigen::VectorXd cgSparse(
const mappedSparseMatrix A,
const mappedVector b
) {
Eigen::ConjugateGradient< mappedSparseMatrix, Eigen::Lower > cg( A ) ;
return cg.solve( b ) ;
}
与R的solve()函数相比:
Comparing with R's solve() function:
B <- matrix( c( 1, 2, 0, 2, 5, 0, 0, 0, 3 ), 3, 3, TRUE )
b <- c( 5, 1, 7 )
B %*% solve( B, b )
[,1]
[1,] 5
[2,] 1
[3,] 7
B %*% cgSparse( as( B, 'dgCMatrix' ), b )
[,1]
[1,] 5
[2,] 1
[3,] 7
这篇关于RcppEigen中的MappedSparseMatrix的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!