Linear models should be estimated using the lm function. In
some cases, lm.fit may be appropriate.
The fastLmPure function provides a reference use case of the GSL
library via the wrapper functions in the RcppGSL package.
The fastLm function provides a more standard implementation of
a linear model fit, offering both a default and a formula interface as
well as print, summary and predict methods.
Lastly, one must be be careful in timing comparisons of
lm and friends versus this approach based on GSL
or Armadillo. The reason that GSL or Armadillo can
do something like lm.fit faster than the functions in
the stats package is because they use the Lapack version
of the QR decomposition while the stats package uses a modified
Linpack version. Hence GSL and Armadillo uses level-3 BLAS code
whereas the stats package uses level-1 BLAS. However,
GSL or Armadillo will choke on rank-deficient model matrices whereas
the functions from the stats package will handle them properly due to
the modified Linpack code. Statisticians want a pivoting scheme of
“pivot only on (apparent) rank deficiency” and numerical
analysts have no idea why statisticians want this so it is not part of
conventional linear algebra software.