Rcpp에서의 회귀 분석에서 대각선 곱하기 RcppEigen

Question

임 RcppEigen을 실행하여 느린 OLS 추정 속도를 높이려고합니다.

나는 Rcpp에서 다음과 같은 olsfunction을 가지고 있지만 아래의 정규 R 경우와 같이 XX 행렬의 대각선을 곱하고 싶습니다. R

ols <- function (y, x) { 
xrd<-crossprod(x) 
xry<-crossprod(x, y) 
diag(xrd)<-1.1*diag(xrd) 
solve(xrd,xry) 

TKS의 P

Answer 1

이 간단 정기

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) 

:

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