[Rcpp-devel] Integrating RcppParallel and RcppEigen

[Joshua Wiley]

Hi JJ,

Thank you, that worked wonderfully.

Thanks to your help, I've been able to add a simple way of bootstrapping
with some success, that seems about 25x faster than direct lm() calls using
4 threads.

I've put it online in case anyone else runs across this and is interested:

https://gist.github.com/JWiley/d9cba55603471f75d438
https://gist.github.com/JWiley/8fe93ae105b1fb244f93

There are some residual issues (the results are not an identical match to
lm(), and occasionally it can make R crash, so I am probably overwriting
something in memory), but once I figure things out, I will fix the gist
files too.

Thanks again for your help and the great package!

Cheers,

Josh







On Tue, Apr 14, 2015 at 7:54 PM, JJ Allaire <jj.allaire at gmail.com> wrote:

> The initialization of your const fields needs to happen within the
> initializer list (rather than inline where the fields are declared).
> If you substitute this code for the code in your example everything
> will compile:
>
>   // QR Decomposition
>   const CPivQR PQR;
>   const Permutation Pmat;
>   const int r;
>
>   // output
>   RMatrix<double> betamat;
>   Eigen::VectorXd betahat;
>
>   // initialize with input and output
>   CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X,
> Rcpp::NumericMatrix betamat)
>     : y(y), X(X), PQR(X), Pmat(PQR.colsPermutation()), r(PQR.rank()),
> betamat(betamat) {}
>
>
>
> On Tue, Apr 14, 2015 at 4:41 AM, Joshua Wiley <jwiley.psych at gmail.com>
> wrote:
> > Hi,
> >
> > I'm interested in combining faster linear model estimation via QR
> > decompositions from RcppEigen with JJs new(er) RcppParallel package to do
> > cross validation and bootstrapping on linear models.
> >
> > As a first step, I've been trying to merge the JSS paper on RcppEigen (
> > http://www.jstatsoft.org/v52/i05/paper ) with about the documentation
> for
> > RcppParallel ( http://rcppcore.github.io/RcppParallel/#examples )
> >
> > ignoring bootstrapping and cross validation for the time being and just
> > getting a linear model to run in the parallel framework.  I've toyed with
> > each separately successfuly, and think I am somewhat close to getting
> them
> > together, but I am running into an error:
> >
> > partest.cpp:28:20: error: 'X' is not a type
> > partest.cpp:29:26: error: 'PQR' is not a type
> > partest.cpp:29:29: error: expected ',' or '...' before '.' token
> > partest.cpp:30:15: error: 'PQR' is not a type
> > partest.cpp:30:18: error: expected ',' or '...' before '.' token
> > partest.cpp: In member function 'virtual void
> CVLM::operator()(std::size_t,
> > std::size_t)':
> > partest.cpp:44:28: warning: comparison between signed and unsigned
> integer
> > expressions [-Wsign-compare]
> > make: *** [partest.o] Error 1
> >
> > My code is below (run via sourceCpp() ), and I can see exactly where the
> > first error occurs, but I don't understand why this is OK when run
> outside
> > of RcppParallel, but so very wrong in conjunction with RcppParallel.
> >
> > Any pointers on either what's wrong or where I should be reading would be
> > deeply appreciated.
> >
> > Josh (code below)
> >
> >
> > // [[Rcpp::depends(RcppParallel)]]
> > // [[Rcpp::depends(RcppEigen)]]
> > #include <Rcpp.h>
> > #include <RcppEigen.h>
> > #include <RcppParallel.h>
> >
> > using namespace std;
> > using namespace Rcpp;
> > using namespace RcppParallel;
> >
> > using Eigen::Lower;
> > using Eigen::Map;
> > using Eigen::MatrixXd;
> > using Eigen::Upper;
> > using Eigen::VectorXd;
> > typedef Map<MatrixXd> MapMatd;
> > typedef Map<VectorXd> MapVecd;
> > typedef Eigen::ColPivHouseholderQR<MatrixXd> CPivQR;
> > typedef CPivQR::PermutationType Permutation;
> >
> > struct CVLM : public Worker
> > {
> >   // design matrix and outcome
> >   const Eigen::VectorXd y;
> >   const Eigen::MatrixXd X;
> >
> >   // QR Decomposition
> >   const CPivQR PQR(X); // ERROR HAPPENS HERE
> >   const Permutation Pmat(PQR.colsPermutation());
> >   const int r(PQR.rank());
> >
> >
> >   // output
> >   RMatrix<double> betamat;
> >   Eigen::VectorXd betahat;
> >
> >   // initialize with input and output
> >   CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X,
> Rcpp::NumericMatrix
> > betamat)
> >     : y(y), X(X), betamat(betamat) {}
> >
> >   void operator()(std::size_t begin, std::size_t end) {
> >
> >     //betahat = PQR.solve(y);
> >     for(int j = begin; j < end; j++) {
> >       betamat(j, j) = X(j, j) + y[j];
> >     }
> >   }
> > };
> >
> > // [[Rcpp::export]]
> > NumericMatrix parallelFit(Eigen::VectorXd dv, Eigen::MatrixXd x) {
> >
> >   // allocate the output matrix
> >   NumericMatrix betamat(x.rows(), x.cols());
> >
> >   // pass input and output
> >   CVLM cvLM(dv, x, betamat);
> >
> >   // parallelFor to do it
> >   parallelFor(0, x.rows(), cvLM);
> >
> >   return betamat;
> > }
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > --
> > Joshua F. Wiley
> > http://joshuawiley.com/
> > ---
> > Postdoctoral Research Fellow
> > Mary MacKillop Institute for Health Research
> > Australian Catholic University
> > ---
> > Senior Analyst, Elkhart Group Ltd.
> > http://elkhartgroup.com
> > Office: 260.673.5518
> >
> > _______________________________________________
> > Rcpp-devel mailing list
> > Rcpp-devel at lists.r-forge.r-project.org
> > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150415/3f93c8c1/attachment.html>