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">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">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><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_extra">this part will always make your code crawl along:<div class="im"><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 class="HOEnZb">
<div class="h5"><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>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">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></blockquote></div><br></div>