[Rcpp-devel] How to increase the coding efficiency

[Douglas Bates]

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/e6b3c442/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: wtls.cpp
Type: text/x-c++src
Size: 494 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20121206/e6b3c442/attachment.cpp>