[Rcpp-devel] How to increase the coding efficiency
wanghonglang2008 at gmail.com
Thu Dec 6 20:43:54 CET 2012
Thanks Douglas, I will study your code and Eigen carefully.
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,
>  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
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the Rcpp-devel