SparseQR for Least Squares

Question

對於我正在構建的應用程序，我需要在大型數據集上運行線性迴歸以獲得殘差。例如，一個數據集的維數超過100萬x 20k。對於較小的數據集，我使用的是RcppArmadillo軟件包中的fastLm，這對於那些目前來說非常適用。隨着時間的推移，這些數據集也將增長超過100萬行。

我的解決方案是使用稀疏矩陣和特徵。我無法找到在RcppEigen中使用SparseQR的好例子。基於很多小時的閱讀（例如：rcpp-gallery,stackoverflow,rcpp-dev mailinglist,eigen docs,rcpp-gallery,stackoverflow以及許多我已經忘記但肯定已經閱讀的內容），我寫了下面這段代碼;

（注：我的第一個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類，它指示支持的減少排序方法。記住

有了這個，我們可以掀起如下：

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