2

the following code yields an error when a Cholesky decomposition is applied to an object of class "matrix"/"array". The problem does not occur when the object is transferred to class "Matrix".

require(Matrix)
set.seed(12345)
nrow <- 200
AM <- matrix(runif(nrow*nrow),nrow,nrow);AM[lower.tri(AM)] <- 0;AM[AM>0.2] <- 0;AM <- AM+t(AM);diag(AM) <- 1000;
A <- Matrix(AM)
cat("class A:",class(A),"\n")
cat("check symmetry A: ",max(abs(A-t(A))),"\n")
cat("rangen eigen A: ",range(eigen(A)$values),"\n")
cat("run chol(A)","\n")
c <- chol(A)
cat("class AM:",class(AM),"\n")
cat("check symmetry AM: ",max(abs(AM-t(AM))),"\n")
cat("rangen eigen AM: ",range(eigen(AM)$values),"\n")
cat("run chol(AM)","\n")
c <- chol(AM)

output:

> Loading required package: Matrix
> > > > > class A: dsCMatrix 
> check symmetry A:  0 
> rangen eigen A:  998.6467 1004.151 
> run chol(A) 
> > class AM: matrix array 
> check symmetry AM:  0 
> rangen eigen AM:  998.6467 1004.151 
> run chol(AM) 
> Error in chol.default(AM) : 
  the leading minor of order 45 is not positive definite
> 

I noticed something odd in other application where solve all of a sudden yielded non-symmetric matrices when factorizing co-variance matrices which were certainly positive -definite. Further, when the factorization was attempted repeatedly the position of the minor changed.

The problem doesn't seem to exist with smaller matrices, say 20x20.

I hope that I overlooked something because otherwise ....

the output of sessioninfo:

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Arch Linux

Matrix products: default
BLAS:   /usr/lib/libopenblasp-r0.3.19.so
LAPACK: /usr/lib/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_AU.UTF-8        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       

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

other attached packages:
[1] Matrix_1.4-0      bit64_4.0.5       bit_4.0.4         data.table_1.14.2

loaded via a namespace (and not attached):
[1] compiler_4.1.2  tools_4.1.2     grid_4.1.2      lattice_0.20-45
> 

Any help much appreciate.

NB:

R version: 4.1.2 (2021-11-01) -- "Bird Hippie"

OS: Linux, Kernel version 5.16.9

CPU: Intel(R) Core(TM) i9-10980XE CPU @ 3.00GHz

Instruction set: avx2, avx512

OPENBLAS_VERSION " OpenBLAS 0.3.19 "

user1407220
  • 353
  • 1
  • 12
  • I can't reproduce the error on macOS or Linux, but I'm also not using OpenBLAS. – Mikael Jagan Feb 19 '22 at 04:06
  • By default, R is built to use internal BLAS/LAPACK libraries. Can you show us the output of `sessionInfo()`? – Mikael Jagan Feb 19 '22 at 04:30
  • Edited ............... But my guess is: it is in openblas. I could reproduce the problem in julia as well, which also relies on openblas. – user1407220 Feb 19 '22 at 05:53
  • You can try building R without an external BLAS to see if it helps. You can refer to the R Installation and Administration [manual](https://cran.r-project.org/doc/manuals/r-devel/R-admin.html) for instructions. – Mikael Jagan Feb 19 '22 at 06:17
  • R's local blas is not a viable option for large linear algebra problems. I have reinstalled R with mkl-blas support ......... that seems to fix it. – user1407220 Feb 19 '22 at 07:19
  • @MikaelJagan According to [this answer](https://stackoverflow.com/a/29986242/6574038), recompiling R is not necessary. – jay.sf Feb 19 '22 at 11:16
  • fyi; both examples run and return the same answer using OpenBLAS 0.3.8 , ubuntu 20.04 – user20650 Feb 19 '22 at 14:13
  • @jay.sf Well, it depends on how the build was configured (and Dirk says the same in that answer). My interpretation of the [BLAS section](https://cran.r-project.org/doc/manuals/r-devel/R-admin.html#BLAS) of the manual is that, if an external BLAS is requested with `--with-blas`, but dynamic linking is _not_ requested with `--enable-blas-shlib`, then BLAS will be statically linked into R. [It seems](https://github.com/archlinux/svntogit-packages/blob/packages/r/trunk/PKGBUILD) that that is precisely how Arch Linux configures its R builds, but I'm not sure. – Mikael Jagan Feb 19 '22 at 16:42
  • @MikaelJagan In Ubuntu 20.04 LTS, doing `sudo apt install libopenblas-base libatlas3-base` to install openblas and atlas, and then `sudo update-alternatives --config libblas.so.3-x86_64-linux-gnu` to select the desired BLAS actually worked fine for me; I switched from blas to openblas and R perfectly started with it. However, OP is using Arch Linux and this probably varies from case to case. BTW, [here a link](https://csantill.github.io/RPerformanceWBLAS/) with benchmarks for different BLAS. mkl appears to run fastest on Intel which is used by OP, whereas openblas fastest for AMD. – jay.sf Feb 20 '22 at 05:01
  • There is a confirmed bug in openblas 0.3.19 related to avx512 – user1407220 Feb 20 '22 at 09:11

1 Answers1

0

Changing blas from openblas to mkl fixed the problem.

there is confirmed bug in openblas 0.3.19 which affects avx512.

user1407220
  • 353
  • 1
  • 12