.Internal(La_rs()) 在某些安装上返回负值，但在其他安装上不返回

Posted 2023-02-18

【中文标题】.Internal(La_rs()) 在某些安装上返回负值，但在其他安装上不返回

【英文标题】：.Internal(La_rs()) returns negative values on some installations but not others

【发布时间】：2020-01-17 22:19:26

【问题描述】：

这是上一个问题的延续： Rfast hd.eigen() returns NAs but base eigen() does not

.Internal(La_rs((x)) 在不同机器上返回不同结果时遇到问题。

我怀疑它可能与数字格式有关，因为在同一台机器上，如果我保存为 CSV 并重新打开，我就不会再收到负数了：

在 Clear Linux 上安装：

> load("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

但是在我的 Windows 安装（以及基于 CentOS 的 HPC 设置）上，我可以直接打开 rdata 文件并且不会得到负值：

> load("input_to_La_rs.Rdata")
> r <- .Internal(La_rs(x, only.values=TRUE))
> sum(r$values < 0)
[1] 0

这与不同的 R 版本/库版本有关吗？一些我不知道的设置？一个错误？

编辑：这是一个更新的示例。它似乎工作不一致，即使在这个特定的安装上，有时我确实得到零：

set.seed(123)
bigm <- matrix(rnorm(2000*2000,mean=0,sd = 3), 2000, 2000)
m <- Rfast::colmeans(bigm)
y <- t(bigm) - m
xx <- crossprod(y)
x <- unname(as.matrix(xx))
b <- .Internal(La_rs(x, TRUE))
sum(b$values < 0)
# [1] 1

另一个更新：事实证明，Rfast 的colmeans 产生的结果与基本 colMeans 略有不同。

    set.seed(123)
    bigm <- matrix(rnorm(2000*2000,mean=0,sd = 3), 2000, 2000)
    m <- colMeans(bigm)
    m <- colmeans(bigm)
    y <- t(bigm) - m
    xx <- crossprod(y)
    x <- unname(as.matrix(xx))
    b <- .Internal(La_rs(x, TRUE))
    sum(b$values < 0)
  # [1] 1

    m <- colMeans(bigm)
    y <- t(bigm) - m
    xx <- crossprod(y)
    x <- unname(as.matrix(xx))
    b <- .Internal(La_rs(x, TRUE))
    sum(b$values < 0)

【问题讨论】：

你能提供一个产生问题的（小）数据集吗？您是从源代码还是从某个存储库安装 R？

我已经在带有 OpenBLAS 的 debian 上尝试了 R 3.6.1，并参考了 BLAS/LAPACK 以及带有 R 提供的 BLAS/LAPACK 的 R devel。在这些环境中，我都没有看到负特征值。但是，sum(eig2$vectors < 0) 对于 OpenBLAS 是 4959，对于参考 BLAS/LAPACK 和 R 的 BLAS/LAPACK 是 4969。也许你可以提供一个Dockerfile 来重现这个问题？

我不再理解这个问题。特征向量中的负分量是预期的！

@RalfStubner 我更新了一个更好的例子，这个确实给出了一个负的特征值。

矩阵x 几乎是秩亏的，即最小的特征值非常接近于零。可能是数值算法变得不稳定并产生了一个小的负值。

【参考方案1】：

Rfast 中的 hd.eigen 函数仅适用于 n

Rfast2 包含一个名为“pca”的函数，适用于任何一种情况，np。也试试那个。在那里，SVD 被有效地执行，从 R 调用“svd”。

【参考方案2】：

我认为差异可以忽略不计。等于或小于 10^(-12) 的矩阵元素的平均差实际上为零。

【讨论】：

它们很小，但它们是下游产生的负特征值和非负特征值之间的差异，最终导致hd.eigen产生NA。 “问题”也是如此，如果有的话，列中的微小差异意味着La_rs 返回负特征值，或者 hd.eigen 将负值转换为 NA 的事实？