Rfast hd.eigen()返回NAs，但base eigen()不返回

Question

我在Rfast中使用hd.eigen时遇到了问题。它给出的结果与包含大多数数据的eigen非常接近，但有时hd.eign会返回空的$vector、NAs或其他不需要的结果。例如：

> 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

除了它破坏了我的脚本之外，这些NAs是否表明我的数据存在更深层次的问题？或者是hd.eig不能处理股票eigen()可以处理的一些情况？一个比另一个好吗？

编辑:根据Ralf的建议，我检查了我的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的高性能计算系统上尝试了相同的例子，但没有得到任何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中生成NAs的表达式为

vectors <- tcrossprod(y, t(FF) * L^(-0.5))

具体地说，L^(-0.5)在索引2000处生成NaN

> L[2000]
[1] -1.136237e-12

然而，在没有返回NAs的两台机器上，L2000是正的(尽管略有不同，在HPC系统上是5.822884e-14，在我的运行Microsoft build of R的Windows机器上是3.022511e-12 )。

编辑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

这能给任何人提供线索吗？这是一个bug吗？

Answer 1

Rfast中的hd.eigen函数仅针对n小于p的情况而设计。