These days I have been working with computation and programming languages. I want to share something with you here.

  1. You cannot expect C++ to magically make your code faster. If speed is of concern, you need profiling to find the bottleneck instead of blind guessing.——Yan Zhou. Thus we have to learn to know how to profile an program in R, Matlab, C++, Python.
  2. When something complicated does not work, I generally try to restart with something simpler, and make sure it works.——Dirk Eddelbuettel.
  3. If you’re calling your function thousands or millions of times, then it might pay to closely examine your memory allocation strategies and figure out what’s temporary.——Christian Gunning.
  4. No, your main issue is not thinking about the computation.  As soon as you write something like
          arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;
    you are in theory land which has very little relationship to practical numerical linear algebra.  If you want to perform linear algebra calculations like weighted least squares you should first take a bit of time to learn about numerical linear algebra as opposed to theoretical linear algebra.  They are very different disciplines.  In theoretical linear algebra you write the solution to a system of linear equations as above, using the inverse of the system matrix.  The first rule of numerical linear algebra is that you never calculate the inverse of a matrix, unless you only plan to do toy examples.  You mentioned sizes of 4000 by 4000 which means that the method you have chosen is doing thousands of times more work than necessary (hint: how do you think that the inverse of a matrix is calculated in practice? – ans: by solving n systems of equations, which you are doing here when you could be solving only one).
    Dirk and I wrote about 7 different methods of solving least squares problems in our vignette on RcppEigen.  None of those methods involve taking the inverse of an n by n matrix.
    R and Rcpp and whatever other programming technologies come along will never be a “special sauce” that takes the place of thinking about what you are trying to do in a computation.——Douglas Bates. |//[[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));
    | }
  5. Repeatedly calling an R function is probably not the smartest thing to do in an otherwise complex and hard to decipher program.—-Dirk Eddelbuettel.
  6. Computers don’t do random things, unlike human beings. Something worked once, is very likely to work whatever times you repeat it as long as the input is the same (unless the function has side effect). So repeating it 1000 times is the same as once.——Yan Zhou
  7. Yan Zhou: Here are a few things people usually do before asking in a mailing list (not just Rcpp list, any such lists like R-help, StackOverflow, etc).
    1. I write a program, it crashes,
    2. I find out the site of crash
    3. I make the program simpler and simpler until it is minimal and the crash is now reproducible.
    4. I still cannot figure out what is wrong with that four or five lines that crash the minimal example
    5. I ask
  8. It does not matter how stupid your questions are. We all asked silly questions before, that is how we learn. But it matters you put in effort to ask the right question. The more effort you put into it, the more specific question you ask and more helpful answers you get.——Yan Zhou.