Karl May
2022-Sep-02 08:29 UTC
[R] stunning matrix multiplication differences Matrix vs. matrix
Hi all, can somebody help me out with an explanation of this:> class(M);class(MM);class(RHS)[1] "dgCMatrix" attr(,"package") [1] "Matrix" [1] "matrix" "array" [1] "dgeMatrix" attr(,"package") [1] "Matrix"> dim(M);dim(MM);dim(RHS)[1] 29418 29418 [1] 29418 29418 [1] 29418 1> max(abs(M-MM))[1] 0> t(RHS)%*%M%*%RHS1 x 1 Matrix of class "dgeMatrix" [,1] [1,] 7264.878> t(RHS)%*%MM%*%RHS1 x 1 Matrix of class "dgeMatrix" [,1] [1,] 5765.366>Apparently the two matrices "M" and "MM" are entirely identical but the difference in results from multiplying both with the same matrix is, say, significant. Even when I transfer M into "matrix" and then back into "Matrix" the differences remain:> M <- as.matrix(M) > t(RHS)%*%M%*%RHS1 x 1 Matrix of class "dgeMatrix" [,1] [1,] 5765.366> M <- Matrix(M) > t(RHS)%*%M%*%RHS1 x 1 Matrix of class "dgeMatrix" [,1] [1,] 7264.878>Frankly, I have no clue what's going on!! This is the session info:> sessionInfo()R version 4.1.3 (2022-03-10) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Arch Linux BLAS/LAPACK: /opt/intel/mkl/lib/intel64/libmkl_gf_lp64.so [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C [3] LC_TIME=C LC_COLLATE=en_AU.UTF-8 [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8 [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C [1] graphics grDevices datasets stats utils parallel methods [8] base [1] pedigreemm_0.3-3 lme4_1.1-28 Matrix_1.4-0 bit64_4.0.5 [5] bit_4.0.4 data.table_1.14.2 loaded via a namespace (and not attached): [1] Rcpp_1.0.8 lattice_0.20-45 MASS_7.3-55 grid_4.1.3 [5] nlme_3.1-155 minqa_1.2.4 nloptr_2.0.0 boot_1.3-28 [9] splines_4.1.3 tools_4.1.3 compiler_4.1.3
Bert Gunter
2022-Sep-03 08:22 UTC
[R] stunning matrix multiplication differences Matrix vs. matrix
FAQ 7.31 ? Your matrices have about 9*10e8 entries after all. I know almost nothing about the Matrix package, but what does crossprod(RHS, Z) %*% RHS give for the two versions of Z = M and MM (assuming Matrix:::crossprod S4 method works here)? Hopefully a knowledgeable Matrix user will correct me if this is wrong. Cheers, Bert On Fri, Sep 2, 2022 at 3:23 PM Karl May <karl.may0 at freenet.de> wrote:> > Hi all, > > can somebody help me out with an explanation of this: > > > class(M);class(MM);class(RHS) > [1] "dgCMatrix" > attr(,"package") > [1] "Matrix" > [1] "matrix" "array" > [1] "dgeMatrix" > attr(,"package") > [1] "Matrix" > > dim(M);dim(MM);dim(RHS) > [1] 29418 29418 > [1] 29418 29418 > [1] 29418 1 > > max(abs(M-MM)) > [1] 0 > > t(RHS)%*%M%*%RHS > 1 x 1 Matrix of class "dgeMatrix" > [,1] > [1,] 7264.878 > > t(RHS)%*%MM%*%RHS > 1 x 1 Matrix of class "dgeMatrix" > [,1] > [1,] 5765.366 > > > > Apparently the two matrices "M" and "MM" are entirely identical but the > difference in results from multiplying both with the same matrix is, say, > significant. > > Even when I transfer M into "matrix" and then back into "Matrix" the > differences remain: > > > M <- as.matrix(M) > > t(RHS)%*%M%*%RHS > 1 x 1 Matrix of class "dgeMatrix" > [,1] > [1,] 5765.366 > > M <- Matrix(M) > > t(RHS)%*%M%*%RHS > 1 x 1 Matrix of class "dgeMatrix" > [,1] > [1,] 7264.878 > > > > Frankly, I have no clue what's going on!! > > This is the session info: > > > sessionInfo() > R version 4.1.3 (2022-03-10) > Platform: x86_64-pc-linux-gnu (64-bit) > Running under: Arch Linux > BLAS/LAPACK: /opt/intel/mkl/lib/intel64/libmkl_gf_lp64.so > > [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C > [3] LC_TIME=C LC_COLLATE=en_AU.UTF-8 > [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8 > [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C > > [1] graphics grDevices datasets stats utils parallel methods > [8] base > [1] pedigreemm_0.3-3 lme4_1.1-28 Matrix_1.4-0 bit64_4.0.5 > [5] bit_4.0.4 data.table_1.14.2 > > loaded via a namespace (and not attached): > [1] Rcpp_1.0.8 lattice_0.20-45 MASS_7.3-55 grid_4.1.3 > [5] nlme_3.1-155 minqa_1.2.4 nloptr_2.0.0 boot_1.3-28 > [9] splines_4.1.3 tools_4.1.3 compiler_4.1.3 > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.