[Rcpp-devel] Detecting Singular Matrices in RcppArmadillo

[Dirk Eddelbuettel]

On 31 January 2022 at 09:13, Zé Vinícius wrote:
| Yes, quoting the paper on which ‘solve’ is based on (
| http://arma.sourceforge.net/armadillo_solver_2020.pdf):
| 
| “The SVD-based solver uses the xGELSD set of functions, which find a
| minimum-norm solution to a linear least squares problem.”

Yes. A longer story though on how this is not what may matter. I should start
by saying that I never found really decent documentation describing this
'split' in the world view as this problem is, in essence, seen differently by
numerical analysis specialists (as for example the LAPACK authors of the
*GELSD functions) and the statistical users which differs in what they focus
on. Which is why R uses a modified version of LINPACK (as I recall going back
to by Ross Ihaka) when computing lm(). And not the LAPACK routines.

Doug was always adamant about this when I wrote the different simple fastLm()
approaches (in RcppGSL, RcppArmadillo, ...) which do _not_ properly account
for rank-deficiency (as R's lm() would) -- Doug also wrote the nicest fastLm
example in RcppEigen.

There is a little bit more in the help pages for the various lmFast()
versions as well as an explicit example (also due to Doug).

Now, I should add that whenever I tried to construct an example on a more
real-world-alike regression problem, I could not come anywhere close to
actually seeing the rank deficiency. But it is unmistakenly there is the
appropriately created case. So buyer beware.

Dirk

-- 
https://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org