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 "