[Rcpp-devel] Detecting Singular Matrices in RcppArmadillo

Very interesting, Dirk. Thanks!

On Mon, Jan 31, 2022 at 12:18 PM Dirk Eddelbuettel <edd at debian.org> wrote:

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

--
Zé Vinícius
https://mirca.github.io