RcppCore / RcppEigen

Rcpp integration for the Eigen templated linear algebra library
Other
110 stars 40 forks source link

Small inconsistency in the `lmBenchmark.R` file? #126

Open statzhero opened 1 year ago

statzhero commented 1 year ago

In the lmBenchmark.R file I read

    ## LDLt Cholesky decomposition with rank detection
    exprs["LDLt"] <- alist(RcppEigen::fastLmPure(mm, y, 2L))

but in the documentation for RcppEigen::fastLm the method is described as

an integer scalar with value [...] 3 for the LDLT Cholesky,

Is this correct?

> sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS 13.3.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] bench_1.1.3         RcppEigen_0.3.3.9.3 speedglm_0.3-5     
[4] biglm_0.9-2.1       DBI_1.1.3           MASS_7.3-58.1      
[7] Matrix_1.5-3       
eddelbuettel commented 1 year ago

It could be a mismatch. The enum that drives is 0-based (in fastLm.cpp)

    enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t, GESDD_t}; 

so that would LLT be 2 and LDLT be 3. The actualy branching is in in the static inline function do_lm() just below. You could add a print statement to be sure... This is all code originally written by @dmbates which has stood some tests of time but it is always possibly that these two got crossed.

statzhero commented 1 year ago

Is it possible that the documentation is behind? It doesn't seem to mention:

    ## SVD using the Lapack subroutine dgesdd and Eigen support
    exprs["GESDD"] <- alist(RcppEigen::fastLmPure(mm, y, 6L))
eddelbuettel commented 1 year ago

[ NNB: You could make this easier by stating which files you take snippets from or use the GitHub feature of a link with an anchor. ]

From what I just quoted, the enum for the GESDD method is six, this line calls with six. Seems correct to me.

statzhero commented 1 year ago

Here is the link to the documentation file if you want to check.

https://github.com/RcppCore/RcppEigen/blob/ab4dcd2df8bde40b26c533896e75d691371fac4a/man/fastLm.Rd#LL35C1-L40C63

\item{method}{an integer scalar with value 0 for the column-pivoted QR decomposition, 1 for the unpivoted QR decomposition, 2 for the LLT Cholesky, 3 for the LDLT Cholesky, 4 for the Jacobi singular value decomposition (SVD) and 5 for a method based on the eigenvalue-eigenvector decomposition of \eqn{\mathbf{X}^\prime\mathbf{X}}{X'X}. Default is zero.}