MappedSparseMatrix in RcppEigen

Question

큰 희소 매트릭스를 해결하기 위해 RcppEigen 패키지로 구현 된 공액 그래디언트 알고리즘을 사용하고 싶습니다.

저는 Rcpp와 C++을 처음 사용하기 때문에 밀도 매트릭스로 시작했습니다.

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

예상대로 작동합니다. 따라서, 나는 이것을 매트릭스를 희소하게하기 위해 확장한다고 생각했다.

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

그러나 이것은 실제로 이상한 결과를 초래합니다. 코드 자체도 오류를주지 않습니다. 누군가가 나를이 오류를 해결하는 데 도움이되기를 바랍니다.

감사합니다.

Answer 1

Eigen :: ConjugateGradient 함수에서 Eigen :: MappedSparseMatrix 형식을 사용해야합니다. R의 해결() 함수와 비교

#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) ; 
} 

:

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