최소 사각형에 대한 SparseQR

Question

응용 프로그램에 대해 나는 잔차를 얻기 위해 큰 데이터 세트에서 선형 회귀 분석을 실행해야합니다. 예를 들어, 하나의 데이터 세트는 차원이 1 백만 x 20k 이상입니다. 더 작은 데이터 세트의 경우 RcppArmadillo 패키지의 fastLm을 사용했습니다. 시간이 지나면 해당 데이터 세트는 1 백만 행을 초과하여 커집니다.

내 솔루션은 희소 행렬과 Eigen을 사용했습니다. RcppEigen에서 SparseQR을 사용하는 좋은 예를 찾을 수 없었습니다. 여러 시간의 독서 (예 : rcpp-gallery, stackoverflow, rcpp-dev mailinglist, eigen docs, rcpp-gallery, stackoverflow 및 그 이상을 잊어 버렸지만 확실하게 읽음)을 기반으로 다음 코드를 작성했습니다.

(NB : 내 첫 번째 C++ 프로그램 - :) 좋은 주시기 바랍니다 - 개선하기 위해 어떤 조언을 환영)

이것은 예를 들어 작동

// [[Rcpp::depends(RcppEigen)]] 
#include <RcppEigen.h> 
using namespace Rcpp; 
using namespace Eigen; 

using Eigen::Map; 
using Eigen::SparseMatrix; 
using Eigen::MappedSparseMatrix; 
using Eigen::VectorXd; 
using Eigen::SimplicialCholesky; 


// [[Rcpp::export]] 
List sparseLm_eigen(const SEXP Xr, 
        const NumericVector yr){ 

    typedef SparseMatrix<double>  sp_mat; 
    typedef MappedSparseMatrix<double> sp_matM; 
    typedef Map<VectorXd>    vecM; 
    typedef SimplicialCholesky<sp_mat> solver; 

    const sp_mat Xt(Rcpp::as<sp_matM>(Xr).adjoint()); 
    const VectorXd Xty(Xt * Rcpp::as<vecM>(yr)); 
    const solver Ch(Xt * Xt.adjoint()); 

    if(Ch.info() != Eigen::Success) return "failed"; 

    return List::create(Named("betahat") = Ch.solve(Xty)); 
} 

;

library(Matrix) 
library(speedglm) 
Rcpp::sourceCpp("sparseLm_eigen.cpp") 

data("data1") 
data1$fat1 <- factor(data1$fat1) 
mm <- model.matrix(formula("y ~ fat1 + x1 + x2"), dat = data1) 

sp_mm <- as(mm, "dgCMatrix") 
y <- data1$y 

res1 <- sparseLm_eigen(sp_mm, y)$betahat 
res2 <- unname(coefficients(lm.fit(mm, y))) 

abs(res1 - res2) 

큰 데이터 세트의 경우 (예상대로) 실패합니다. 나의 초기 의도는 SparseQR을 해결사로 사용하는 것이지만이를 구현하는 방법을 모르겠습니다.

그래서 내 질문 - 누군가 RcppEigen과 스파 스 매트릭스에 대한 QR 분해를 구현하는 데 도움이 될 수 있습니까?

Answer 1

Eigen으로 희소 솔버를 쓰는 방법은 조금 일반적입니다. 이것은 주로 스파 스 솔버 클래스가 훌륭하게 설계 되었기 때문입니다. 그들은 guide explaining their sparse solver classes를 제공합니다. 질문은 SparseQR에 초점을 맞추기 때문에 문서에서는 Solver를 초기화하는 데 필요한 두 개의 매개 변수 인 SparseMatrix 클래스 유형과 지원되는 채우기 감소 순서 지정 방법을 나타내는 OrderingMethods 클래스가 있음을 나타냅니다. 이와

염두에두고, 우리는 다음과 채찍 수 있습니다

// [[Rcpp::depends(RcppEigen)]] 
#include <RcppEigen.h> 
#include <Eigen/SparseQR> 

// [[Rcpp::export]] 
Rcpp::List sparseLm_eigen(const Eigen::MappedSparseMatrix<double> A, 
          const Eigen::Map<Eigen::VectorXd> b){ 

    Eigen::SparseQR <Eigen::MappedSparseMatrix<double>, Eigen::COLAMDOrdering<int> > solver; 
    solver.compute(A); 
    if(solver.info() != Eigen::Success) { 
    // decomposition failed 
    return Rcpp::List::create(Rcpp::Named("status") = false); 
    } 
    Eigen::VectorXd x = solver.solve(b); 
    if(solver.info() != Eigen::Success) { 
    // solving failed 
    return Rcpp::List::create(Rcpp::Named("status") = false); 
    } 

    return Rcpp::List::create(Rcpp::Named("status") = true, 
          Rcpp::Named("betahat") = x); 
} 

참고 : 여기에 우리가 항상 처음를 확인해야 명명 된 status 변수를 전달하는 목록을 만들 수 있습니다. 이는 수렴이 분해와 해결의 두 영역에서 발생하는지 여부를 나타냅니다. 모두 체크 아웃하면 betahat 계수를 전달합니다.

---

테스트 스크립트 :

library(Matrix) 
library(speedglm) 
Rcpp::sourceCpp("sparseLm_eigen.cpp") 

data("data1") 
data1$fat1 <- factor(data1$fat1) 
mm <- model.matrix(formula("y ~ fat1 + x1 + x2"), dat = data1) 

sp_mm <- as(mm, "dgCMatrix") 
y <- data1$y 

res1 <- sparseLm_eigen(sp_mm, y) 
if(res1$status != TRUE){ 
    stop("convergence issue") 
} 
res1_coef = res1$betahat 
res2_coef <- unname(coefficients(lm.fit(mm, y))) 

cbind(res1_coef, res2_coef) 

출력 :

 res1_coef res2_coef 
[1,] 1.027742926 1.027742926 
[2,] 0.142334262 0.142334262 
[3,] 0.044327457 0.044327457 
[4,] 0.338274783 0.338274783 
[5,] -0.001740012 -0.001740012 
[6,] 0.046558506 0.046558506