c++ - 最小二乘的 SparseQR

Question

对于我正在构建的应用程序，我需要对大型数据集运行线性回归以获得残差。例如，一个数据集的维度超过 100 万 x 20k。对于较小的数据集，我使用的是 RcppArmadillo 包中的 fastLm - 目前非常适合这些数据集。随着时间的推移，这些数据集也将增长到超过 100 万行。

我的解决方案是使用稀疏矩阵和 Eigen。我找不到在 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 分解吗？

---

我使用的解决方案

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

根据@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);
}

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