在Rcpp中将ols回归中的对角相乘
在Rcpp中将ols回归中的对角相乘
提问于 2014-04-23 11:38:23
我试图通过运行RcppEigen来加速我缓慢的OLS估计。
我在Rcpp中有以下的函数，但是想把XX矩阵中的对角线相乘，就像下面的正则R情况一样？
fastolsCpp <- '
const MapMatd X(as<MapMatd>(XX));
const MapVecd y(as<MapVecd>(yy));
const LLT<MatrixXd> llt(AtA(X));
const VectorXd betahat(llt.solve(X.adjoint() * y));
return wrap(betahat);
'
fastols <- cxxfunction(signature(XX = "matrix", yy = "numeric"), fastolsCpp, "RcppEigen",incl)正则R
ols <- function (y, x) {
xrd<-crossprod(x)
xry<-crossprod(x, y)
diag(xrd)<-1.1*diag(xrd)
solve(xrd,xry)Tks P
回答 1
Stack Overflow用户
回答已采纳
发布于 2014-04-23 17:55:38
这很简单：
library(RcppEigen)
library(inline)
set.seed(42)
x <- 1:10
y <- 3+2*x+rnorm(10)
incl <- '
using Eigen::LLT;
using Eigen::Lower;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;
using Rcpp::as;
typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
typedef Eigen::Map<Eigen::VectorXd> MapVecd;
inline MatrixXd AtA(const MapMatd& A) {
int n(A.cols());
return MatrixXd(n,n).setZero().selfadjointView<Lower>().rankUpdate(A.adjoint());
}
'
fastolsCpp <- '
const MapMatd X(as<MapMatd>(XX));
const MapVecd y(as<MapVecd>(yy));
MatrixXd A=AtA(X);
const int k=A.rows();
for (int i=0; i<k; i++) {
A(i,i) *= 1.1;
}
const LLT<MatrixXd> llt(A);
const VectorXd betahat(llt.solve(X.adjoint() * y));
return wrap(betahat);
'
fastols <- cxxfunction(signature(XX = "matrix", yy = "numeric"), fastolsCpp, "RcppEigen", incl)
all.equal(fastols(cbind(1,x), y),
c(ols(y, cbind(1,x))))
#[1] TRUE
library(microbenchmark)
microbenchmark(ols(y, cbind(1,x)), fastols(cbind(1,x), y))
# Unit: microseconds
# expr min lq median uq max neval
# ols(y, cbind(1, x)) 137.585 139.1775 140.571 148.2945 359.642 100
# fastols(cbind(1, x), y) 12.533 13.4830 15.904 16.1720 55.743 100页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/23243166
复制相关文章
相似问题
腾讯云开发者
