对于我正在构建的应用程序，我需要在大型数据集上运行线性回归以获得残差。例如，一个数据集的维度超过100万x 20 k。对于较小的数据集，我使用RcppArmadillo包中的fastLm-目前对这些数据集效果很好。随着时间的推移，这些数据集也将超过100万行。
我的解决方案是使用稀疏矩阵和特征值。我无法找到一个在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分解吗？

我使用的溶液

• 免责声明：* 我不知道它是否正确-它适用于我的特定问题，似乎提供了正确的结果。

根据@TomWenseleers的评论添加我的解，这里是。我不能让Eigen's LeastSquareConjugateGradient工作。我想是因为文档声明A'A必须是正定的。所以我没有求解A = bx，而是使用Eigen's Conjugate Gradient和对角预处理器求解A'A = A'bx。然后使用系数计算拟合值和残差。

#include <RcppEigen.h>
#include <Eigen/IterativeLinearSolvers>
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Rcpp::List sparse_cg(const Eigen::Map<Eigen::SparseMatrix<double> > &A,
                     const Eigen::Map<Eigen::VectorXd> &b,
                     const Eigen::Map<Eigen::VectorXd> &x0,
                     const int &maxiter,
                     const double &tol) {
    
  // Set-up the Conjugate Gradient solver
  Eigen::ConjugateGradient<Eigen::SparseMatrix<double>,
                           Eigen::Lower|Eigen::Upper,
                           Eigen::DiagonalPreconditioner<double> > cg;
  // Initialize solver
  cg.setMaxIterations(maxiter);
  cg.setTolerance(tol);
  cg.compute(A);

  // Solve system with guess (i.e. old solutions)
  Eigen::VectorXd coef = cg.solveWithGuess(b, x0);

  // Solver convergence stats
  int iter = cg.iterations();
  double err = cg.error();

  return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
                            Rcpp::Named("itr") = iter,
                            Rcpp::Named("error") = err);
}