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.tri...@gmail.com>wrote:
> this part will always make your code crawl along:
>
>
> On Wed, Dec 5, 2012 at 11:09 AM, Honglang Wang <wanghonglang2...@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@lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
typedef Eigen::MatrixXd Mat;
typedef Eigen::Map<Mat> MMat;
typedef Eigen::HouseholderQR<Mat> QR;
typedef Eigen::VectorXd Vec;
typedef Eigen::Map<Vec> MVec;
// [[Rcpp::export]]
Rcpp::List wtls(const MMat X, const MVec y, const MVec sqrtwts) {
return Rcpp::List::create(Rcpp::Named("betahat") =
QR(sqrtwts.asDiagonal()*X).solve(sqrtwts.asDiagonal()*y));
}
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
