SparseQR for Least Squares

Question

私が構築しているアプリケーションでは、残差を得るために大きなデータセットに対して線形回帰を実行する必要があります。たとえば、1つのデータセットは100万x 20,000次元以上です。小さなデータセットの場合は、RcppArmadilloパッケージのfastLmを使用していました。時間が経つにつれて、それらのデータセットも100万行を超えるようになります。

私の解決策は、スパース行列と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に焦点を当てているので、ソルバーを初期化するために必要なパラメータはSparseMatrixクラスタイプとサポートされているフィルリダクションの順序付け方法を指示するOrderingMethodsクラスの2つです。これにより

心の中で、私たちは次のことをかき立てることができます。

// [[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変数を渡すリストを作成します。これは、収束が分解と解決の2つの領域で起こるかどうかを示します。すべてチェックアウトすると、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