[Rcpp-devel] Integrating RcppParallel and RcppEigen

JJ Allaire jj.allaire at gmail.com
Tue Apr 14 23:57:02 CEST 2015


Yes, we'd LOVE to see this turned into an Rcpp Gallery article once
you've got everything ironed out!

Details on creating Gallery articles here:
https://github.com/jjallaire/rcpp-gallery/wiki/Contributing-to-the-Rcpp-Gallery

On Tue, Apr 14, 2015 at 5:36 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
> 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
>
>


More information about the Rcpp-devel mailing list