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

Question

我无法在R中找到任何函数或程序包来计算bigmemory（来自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

- Mahin

---

这些中有哪个可用 ?qr，或者 ?Matrix::qr，或者 ?MASS::Null？ - user20650

是的，我知道这些函数，但它们不能用于大矩阵（S4类），或者我无法将它们用于大矩阵。我只能将这些函数用于普通矩阵，而不能用于大矩阵。 - Mahin

bigalgebra 包已经开始了一些方法，但 QR 函数是不完整的。然而，这个分支 https://github.com/cdeterman/bigalgebra 添加了 QR 功能。它会发出警告 警告：这是不建议的 - 你可以问作者为什么。 - user20650

我有一个大矩阵。例如：library(bigmemory) a<-big.matrix(1000000, 1000, type='double', init=5) options(bigmemory.allow.dimnames=TRUE) 但是我找不到它的零空间或QR分解。谢谢。 - Mahin

@F. Privé 是的，我真的想对一个big.matrix进行完整的QR分解（或完全SVD），然后进入它的零空间。 - Mahin

Answer 1

如果您想计算矩阵的完整SVD，您可以使用软件包bigstatsr按块执行计算。 FBM代表文件支持的大型矩阵，类似于软件包bigmemory的文件支持的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分钟。

- F. Privé

---

你好。谢谢你的回答，但我需要计算一个 m×n 矩阵 A 的“完全”奇异值分解（A=u.d.t(v)），其中 u 是一个 m×m 矩阵，d 是一个 m×n 矩阵，v 是一个 n×n 矩阵。我真的需要矩阵 v 才能进入零空间，但是从你提出的方法得到的矩阵元素与我想要的不同。例如，请参见：https://math.stackexchange.com/questions/1771013/null-space-from-svd - Mahin

@Mahin 从我的回答中，你可以得到 u 是 m x n，d 是对角线 n x n（这里只给出了对角线），而且 v 就像你想要的那样是 n x n。从这三个矩阵完全可以完美地重构出 a。 - F. Privé

据我所理解，你提到的完整奇异值分解只是在u中添加了一些（无用）列。v应该是相同的，您可以通过比较具有更多行而不是列的任何矩阵的svd(mat)$v和svd(mat, nu = nrow(mat), nv = ncol(mat))$v来验证它。 - F. Privé

@ F. Privé 您是正确的。我错了。在这两种情况下，svd（mat）$v和svd（mat，nu = nrow（mat），nv = ncol（mat））$v，矩阵元素是相同的，现在我可以从v计算零空间。非常感谢您。 - Mahin

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

@Mahin 在这里，这个程序可以工作是因为你有一个相对较小的维度（1000）。 - F. Privé

Answer 2

@Mahon @user20650 @F.Privė 为了更清楚地表达，我联系了BigMemory团队并询问：

基本上，是否有一个能够处理大内存矩阵的QR函数(QR分解)实现？

我觉得搞清楚最初的问题很有用。@F.Privė - 很好的回答。希望你的回答和他们的回复能够帮助指导未来的人们。他们的回复如下：

感谢您的留言。目前还没有QR分解的实现。理想情况下，如果矩阵是密集的，您可以使用Householder反射；如果是稀疏的，则可以使用Givens旋转。

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

请注意，如果您的矩阵是(高且瘦或短而胖)，那么SO解法是可以的。然而，每次求解叉积时都会失去一些数值稳定性。如果您计划反转矩阵，这可能是个问题。

- Technophobe01

---

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

感谢Techno的有用评论。所以答案真的是没有QR分解或空间计算的开箱即用，因此需要编写C++代码。 - user20650

@user20650 - 是的，那是目前的情况 - 很不幸。Michael K.在这里非常有帮助（又名BigMemory）。话虽如此，我确实喜欢F.Privė的替代方法，我计划更详细地研究一下这个包。 - Technophobe01

@ techno 感谢您的关注。我目前正在使用F.Privė的解决方案，但我希望相关功能很快就能被编写出来。 - Mahin