在 R 中计算大矩阵的零空间

Posted 2023-03-11

【中文标题】在 R 中计算大矩阵的零空间

【英文标题】：Computing the null space of a bigmatrix in R

【发布时间】：2018-02-25 11:17:27

【问题描述】：

我找不到任何函数或包来计算 R 中 bigmatrix（来自 library(bigmemory)）的空空间或（QR 分解）。例如：

library(bigmemory)

a <- big.matrix(1000000, 1000, type='double', init=0)

我尝试了以下但得到了显示的错误。如何找到bigmemory 对象的空空间？

a.qr <- Matrix::qr(a)
# Error in as.vector(data) : 
#   no method for coercing this S4 class to a vector
q.null <- MASS::Null(a)
# Error in as.vector(data) : 
#   no method for coercing this S4 class to a vector

【问题讨论】：

做这些工作?qr，或?Matrix::qr，或?MASS::Null

是的。我有，但是这些函数不适用于大矩阵（S4 类），或者我不能将它们用于大矩阵。我只能将这些函数用于常规矩阵，不能用于大矩阵。

好吧，我不确定你是否有一个大的 矩阵或 bigmatrix ;)。目前，您的问题是off topic，因为它直接要求包推荐，并且在目前的状态下它可能会被关闭。但这很有趣。请您编辑您的问题并提供更多详细信息。例如，您能否添加一个 bigmatrix 的小示例（包括使用的任何包），说明标准工具如何不起作用，并可能要求提供替代方案。谢谢

bigalgebra 包已经启动了一些方法，但是 QR 功能不完整，但是这个 fork ，github.com/cdeterman/bigalgebra 增加了 QR 功能。它确实给出了警告警告：不建议这样做。 - 你可以问作者为什么

我有一个大矩阵。例如： library(bigmemory) a

【参考方案1】：

如果要计算矩阵的完整 SVD，可以使用包 bigstatsr 进行分块计算。 FBM 代表 Filebacked Big Matrix，是一个类似于 bigmemory 包的 filebacked big.matrix 对象的对象。

library(bigstatsr)
options(bigstatsr.block.sizeGB = 0.5)

# Initialize FBM with random numbers
a <- FBM(1e6, 1e3)
big_apply(a, a.FUN = function(X, ind) 
  X[, ind] <- rnorm(nrow(X) * length(ind))
  NULL
, a.combine = 'c')

# Compute t(a) * a
K <- big_crossprodSelf(a, big_scale(center = FALSE, scale = FALSE))

# Get v and d where a = u * d * t(v) the SVD of a
eig <- eigen(K[])
v <- eig$vectors
d <- sqrt(eig$values)

# Get u if you need it. It will be of the same size of u
# so that I store it as a FBM.
u <- FBM(nrow(a), ncol(a))
big_apply(u, a.FUN = function(X, ind, a, v, d) 
  X[ind, ] <- sweep(a[ind, ] %*% v, 2, d, "/")
  NULL
, a.combine = 'c', block.size = 50e3, ind = rows_along(u),
a = a, v = v, d = d)

# Verification
ind <- sample(nrow(a), 1000)
all.equal(a[ind, ], tcrossprod(sweep(u[ind, ], 2, d, "*"), v))

这在我的电脑上大约需要 10 分钟。

【讨论】：

嗨。感谢您的回答，但我需要计算 m×n 矩阵 A (A=udt(v)) 的“完整”SVD​​，其中 u 是 m×m 矩阵，d 是 m×n 矩阵，v 是一个 n×n 矩阵。我真的需要矩阵 v 进入零空间，但是从您提出的方法中获得的矩阵元素与我想要的不同。例如参见：math.stackexchange.com/questions/1771013/null-space-from-svd

@Mahin 从我的回答中，你会得到 u 是 m x n，d 是对角线 n x n（这里只给出对角线），v 是你想要的 n x n。从这三个矩阵可以完美重构a的意义上来说，它是完整的。

@Mahin 据我了解，您所说的完整 SVD 只是在u 中添加一些（无用的）列。 v 应该相同，您可以通过比较 svd(mat)$v 和 svd(mat, nu = nrow(mat), nv = ncol(mat))$v 来验证行数多于列数的任何矩阵。

@F. Privé 你是对的。我错了。在这两种情况下，svd(mat)$v 和 svd(mat, nu = nrow(mat), nv = ncol(mat))$v，矩阵元素是相同的，我现在可以从 v 计算空空间。谢谢你非常喜欢。

是否有可能在您的计算机上拥有一个 1e6x1e6 大矩阵并通过 R 中的这些命令计算其 SVD？我不知道如何处理这种大小的矩阵。谢谢。

【参考方案2】：

@Mahon @user20650 @F.Privė 为了清楚起见，我联系了 bigmemory 团队并询问了

本质上，是否存在适用于大内存矩阵的 QR 函数（QR 分解）的实现？

我觉得弄清楚最初提出的问题很有用。 @F.Privė - 很好的答案。希望您的回答和他们的回答将有助于将来指导人们。他们的回应如下：

感谢您的来信。目前没有 qr 分解的实现。理想情况下，您可以使用 Householder 反射（如果矩阵密集）或 Givens 旋转（如果矩阵稀疏）来实现这一点。

irlba 包 与 bigmemory 兼容。它提供了截断的奇异值分解。因此，如果您的矩阵相对稀疏，您可以在矩阵的秩处截断。这可能是您最好的选择。如果您不知道排名，那么您可以使用包迭代更新截断。

请注意，如果您的矩阵是（又高又瘦或又矮又胖），那么 SO 解决方案是可以的。然而，任何时候你求助于计算叉积，你都会失去一些数值稳定性。如果您打算反转矩阵，这可能是一个问题。

【讨论】：

你好。非常感谢您的指导。但我不想进行线性模型拟合。我只想对 big.matrix 进行完整的 QR 分解（或完整的 SVD），然后进入它的空空间。我运行了“RcppEigen”包，但这个包没有给我 mat.obj 的 QR 分解，只给出了线性模型拟合。谢谢。

感谢 techno，有用的 cmets。所以答案真的是没有 QR 分解，或开箱即用的空空间计算，因此需要编写 c++ 代码。

@user20650 - 是的，这就是目前的情况 - 不幸的是。 Michael K. 在这里非常有帮助（A.k.a BigMemory）。也就是说，我确实喜欢 F.Privė 替代方法，我计划更详细地了解一个包。

@techno 感谢您的关注。我目前正在做F.Privė的方案，但希望相关功能快点写出来。