Strange problem with cholesky in matrix class .............. possible bug?

Question

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 "

Answer 1

Changing blas from openblas to mkl fixed the problem.

there is confirmed bug in openblas 0.3.19 which affects avx512.