If you're open to using Rcpp and its related packages, the existing examples for the fastLm()
show you how to do this with
- Eigen via RcppEigen
- Armadillo via RcppArmadillo
- the GSL via RcppGSL
where the latter two will use the same BLAS as R, and Eigen has something internal that can often be faster. All the packages implement lm()
using the decompositions offered by the language (often using solve or related, but switching to SVD is straightforward once you have the toolchain working for you).
Edit: Here is an explicit example. Use the following C++ code:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::vec getEigen(arma::mat M) {
return arma::eig_sym(M);
}
save in a file "eigenEx.cpp". Then all it takes is this R code:
library(Rcpp) ## recent version for sourceCpp()
sourceCpp("eigenEx.cpp") ## converts source file into getEigen() we can call
so that we can run this:
set.seed(42)
X <- matrix(rnorm(3*3), 3, 3)
Z <- X %*% t(X)
getEigen(Z) ## calls function created above
and I get the exact same eigen vector as from R. It really doesn't get much easier.
It also lets Armadillo chose what method to use for the Eigen decomposition, and as David hinted, this is something quicker than a full-blown SVD (see the Armadillo docs).