我在使用hd.eigen
in时遇到问题Rfast
。它给出了与大多数数据非常接近的结果eigen
,但有时 hd.eign 会返回一个空的 NAs$vector
或其他不需要的结果。例如:
> set.seed(123)
> bigm <- matrix(rnorm(2000*2000,mean=0,sd = 3), 2000, 2000)
>
> e3 = eigen(bigm)
> length(e3$values)
[1] 2000
> length(e3$vectors)
[1] 4000000
> sum(is.na(e3$vectors) == TRUE)
[1] 0
> sum(is.na(e3$vectors) == FALSE)
[1] 4000000
>
> e4 = hd.eigen(bigm, vectors = TRUE)
> length(e4$values)
[1] 2000
> length(e4$vectors)
[1] 4000000
> sum(is.na(e4$vectors) == TRUE)
[1] 2000
> sum(is.na(e4$vectors) == FALSE)
[1] 3998000
除了它破坏了我的脚本之外,这些 NA 是否表明我的数据存在更深层次的问题?还是hd.eig
无法处理股票eigen()
可以处理的某些情况?这个比那个好吗?
编辑:根据拉尔夫的建议,我检查了我的 BLAS 版本,看起来 R 可能正在寻找错误的版本/在错误的位置:
~ $ ldd /usr/lib64/R/bin/exec/R
linux-vdso.so.1 (0x00007ffeec3b9000)
libR.so => not found
libRblas.so => not found
libgomp.so.1 => /usr/lib64/libgomp.so.1 (0x00007feb27ef2000)
libpthread.so.0 => /usr/lib64/libpthread.so.0 (0x00007feb27ecf000)
libc.so.6 => /usr/lib64/libc.so.6 (0x00007feb27cdb000)
/lib64/ld-linux-x86-64.so.2 => /usr/lib64/ld-linux-x86-64.so.2 (0x00007feb27f7b000)
另外,我不清楚 openBLAS 是否等同于其他发行版中默认安装的 BLAS。
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-generic-linux-gnu (64-bit)
Running under: Clear Linux OS
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblas_nehalemp-r0.3.6.so
编辑 2:我在基于 CentOS 的 HPC 系统上尝试了相同的示例,但没有得到任何 NA。在那里,sessionInfo()
揭示了:
BLAS/LAPACK: /hpc/packages/minerva-centos7/intel/parallel_studio_xe_2019/compilers_and_libraries_2019.0.117/linux/mkl/lib/intel64_lin/libmkl_gf_lp64.so
编辑3:hd.eign
产生NA的表达式是
vectors <- tcrossprod(y, t(FF) * L^(-0.5))
具体来说,L^(-0.5)
在索引 2000 处产生 NaN
> L[2000]
[1] -1.136237e-12
但是,在没有返回 NA 的两台机器上,L[2000] 是正数(虽然略有不同,但5.822884e-14
在 HPC 系统和3.022511e-12
我的运行 Microsoft 构建的 R.
编辑 4:差异似乎源于基函数,它从问题机器上的矩阵eigen()
返回一个负值,而不是其他两个。我保存了对象并在计算机之间打开,所以我知道输入是完全相同的。crossprod()
xx
xx
eigen()
编辑5:我钻了一层更深,发现原来的负值来自这个语句eigen()
z <- if (!complex.x)
.Internal(La_rs(x, only.values))
else .Internal(La_rs_cmplx(x, only.values))
编辑 6:如果我保存为 CSV 然后重新打开,问题计算机不会产生负特征值。
> load("/home/james/nfs-cloud/PanosLab/CircRNA/input_to_La_rs.Rdata")
> r <- .Internal(La_rs(as.matrix(x), only.values = FALSE))
> sum(r$values < 0)
[1] 1
> write.csv(x, "test_for_internal.csv", row.names = FALSE)
> x <- read.csv("test_for_internal.csv")
> r <- .Internal(La_rs(as.matrix(x), only.values = FALSE))
> sum(r$values < 0)
[1] 0
这会给任何人一个线索吗?这是一个错误吗?