[Rcpp-devel] How to increase the coding efficiency

[Honglang Wang]

Thanks Douglas, I will study your code and Eigen carefully.

Best wishes!

Honglang Wang

Office C402 Wells Hall
Department of Statistics and Probability
Michigan State University
1579 I Spartan Village, East Lansing, MI 48823
wangho16 at msu.edu



On Thu, Dec 6, 2012 at 1:33 PM, Douglas Bates <bates at stat.wisc.edu> wrote:

> I don't know as much about Armadillo as I do about Eigen so I will cheat
> and write using RcppEigen instead of RcppArmadillo.  Page 6 of the Eigen
> tutorial at http://eigen.tuxfamily.org/dox/ discusses decompositions and
> solving linear systems of equations.  One of the simplest ways of solving a
> weighted least squares system, if you only want the solution itself, is
> with a QR decomposition of which there are 3 in Eigen.  The simplest is
> HouseholderQR (named after A.S. Householder who first described an
> elementary reflector).  For a least squares solution (assuming the model
> matrix has full column rank) you simply use the solve method. For a
> weighted least squares solution you premultiply both the model matrix and
> the response by the square roots of the weights.
>
> It can be collapsed to a one-liner, as in the enclosed.  A test that it
> works for unweighted least squares (i.e. using unit weights) is
>
> > sourceCpp("/tmp/wtls.cpp")
> > coef(lm(optden ~ carb, Formaldehyde))(Intercept)        carb
> 0.005085714 0.876285714
> > with(Formaldehyde, wtls(model.matrix(~carb), optden, rep.int(1,
> length(optden))))
> $betahat
> [1] 0.005085714 0.876285714
>
> The use of the other decompositions (Cholesky, ColPivHouseholderQR,
> Symmetric Eigendecomposition, SVD) are described in the RcppEigen vignette.
>
> By the way, the sourceCpp, evalCpp and cppFunction capabilities developed
> by the Rcpp authors is very close to magic.  They are to be congratulated
> on a remarkable accomplishment.
>
>
> On Wed, Dec 5, 2012 at 1:16 PM, Tim Triche, Jr. <tim.triche at gmail.com>wrote:
>
>> this part will always make your code crawl along:
>>
>>
>> On Wed, Dec 5, 2012 at 11:09 AM, Honglang Wang <
>> wanghonglang2008 at gmail.com> wrote:
>>
>>>   arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;
>>>
>>
>> first time I wrote a GLM engine, I wrote it the way statistics books
>> illustrate it (i.e. actually invert things to do IWLS) and it crawled.  I
>> took it to a physicist friend who went through alternate stages of disgust
>> and laughter then told me never to invert anything I didn't have to.
>>
>> You should take Doug Bates' advice, it could save you a great deal of
>> time.  Never bit fiddle when you can use a better algorithm.
>>
>>
>> --
>> *A model is a lie that helps you see the truth.*
>> *
>> *
>> Howard Skipper<http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf>
>>
>>
>> _______________________________________________
>> 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/20121206/78258caa/attachment-0001.html>