Thanks Douglas, I will study your code and Eigen carefully.<br><br clear="all"><div>Best wishes!</div><div> </div><div>Honglang Wang</div><div> </div><div>Office C402 Wells Hall</div><div>Department of Statistics and Probability</div>
<div>Michigan State University</div><div>1579 I Spartan Village, East Lansing, MI 48823</div><div><a href="mailto:wangho16@msu.edu" target="_blank">wangho16@msu.edu</a></div><br>
<br><br><div class="gmail_quote">On Thu, Dec 6, 2012 at 1:33 PM, Douglas Bates <span dir="ltr"><<a href="mailto:bates@stat.wisc.edu" target="_blank">bates@stat.wisc.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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 <a href="http://eigen.tuxfamily.org/dox/" target="_blank">http://eigen.tuxfamily.org/dox/</a> 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.<div>
<br></div><div>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</div><div><br></div><div><div>> sourceCpp("/tmp/wtls.cpp")</div>
<div>> coef(lm(optden ~ carb, Formaldehyde))(Intercept) carb <br></div><div>0.005085714 0.876285714 </div><div>> with(Formaldehyde, wtls(model.matrix(~carb), optden, <a href="http://rep.int" target="_blank">rep.int</a>(1, length(optden))))</div>
<div>$betahat</div><div>[1] 0.005085714 0.876285714</div></div><div><br></div><div>The use of the other decompositions (Cholesky, ColPivHouseholderQR, Symmetric Eigendecomposition, SVD) are described in the RcppEigen vignette.</div>
<div><br></div><div>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.</div><div class="gmail_extra">
<br><br><div class="gmail_quote"><div><div class="h5">On Wed, Dec 5, 2012 at 1:16 PM, Tim Triche, Jr. <span dir="ltr"><<a href="mailto:tim.triche@gmail.com" target="_blank">tim.triche@gmail.com</a>></span> wrote:<br>
</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5">
<div class="gmail_extra">this part will always make your code crawl along:<div><br><br><div class="gmail_quote">On Wed, Dec 5, 2012 at 11:09 AM, Honglang Wang <span dir="ltr"><<a href="mailto:wanghonglang2008@gmail.com" target="_blank">wanghonglang2008@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div> arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;<br></div><div></div></blockquote></div>
<div class="gmail_extra"><br></div></div>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. </div>
<div class="gmail_extra"><br></div><div class="gmail_extra">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.</div><div>
<div><div class="gmail_extra">
<br></div><div class="gmail_extra"><br>-- <br><i>A model is a lie that helps you see the truth.</i><div><i><br></i><div><a href="http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf" target="_blank">Howard Skipper</a></div>
</div><br>
</div>
</div></div><br></div></div><div class="im">_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></div></blockquote></div><br></div>
</blockquote></div><br>