对于我正在构建的应用程序,我需要对大型数据集运行线性回归以获得残差。例如,一个数据集的尺寸超过100万x 20k。对于较小的数据集,我使用的是fastLm
包中的RcppArmadillo
-目前非常适合这些数据集。随着时间的推移,这些数据集还将增长到超过一百万行。
我的解决方案是使用稀疏矩阵和本征。我找不到在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分解吗?
最佳答案
如何使用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
变量,该变量应首先检查。这表明收敛是否发生在两个 Realm :分解和求解。如果全部 checkout ,则我们传递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