Multiplying the diagonal in an ols regression in Rcpp RcppEigen

Question

Im trying to speed up my slow OLS estimation by running RcppEigen.

I have the following olsfunction in Rcpp but would like to multiply the diagonal in the XX matrix like in the regular R case below?

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)

Regular 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

This is straightforward:

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