大型稀疏矩阵的 Moore-Penrose 广义逆

Posted 2023-03-11

【中文标题】大型稀疏矩阵的 Moore-Penrose 广义逆

【英文标题】：Moore-Penrose generalized inverse of a large sparse matrix

【发布时间】：2014-07-12 14:15:14

【问题描述】：

我有一个有几万行和几万列的方阵，除了大量的0 之外只有几个1，所以我使用Matrix 包以有效的方式将其存储在R 中。 base::matrix 对象无法处理这么多的单元格，因为我的内存不足。

我的问题是我需要此类矩阵的 逆 以及 Moore-Penrose 广义逆，目前我无法计算。

我尝试过的：

solve 产生 Error in LU.dgC(a) : cs_lu(A) failed: near-singular A (or out of memory) 错误 MASS::ginv 与 Matrix 类不兼容 没有直接的方法可以将稀疏的Matrix 转换为例如bigmemory::big.matrix，后者无论如何都不适用于MASS::ginv

如果我尝试计算矩阵的 Choleski 分解以稍后调用 Matrix::chol2inv，我会收到以下错误消息：

Error in .local(x, ...) : 
  internal_chm_factor: Cholesky factorization failed
In addition: Warning message:
In .local(x, ...) :
  Cholmod warning 'not positive definite' at file ../Cholesky/t_cholmod_rowfac.c, line 431

基于related question，我还尝试了单个节点上的pbdDMAT 包，但pbdDMAT::chol 产生了Cholmod error 'out of memory' at file ../Core/cholmod_memory.c, line 147 错误消息

简而言之问题：除了在具有大量 RAM 的计算机上使用matrix 类之外，还有什么方法可以计算这种稀疏矩阵的逆矩阵和 Moore-Penrose 广义逆矩阵?

---

快速可重现的示例：

library(Matrix)
n <- 1e5
m <- sparseMatrix(1:n, 1:n, x = 1)
m <- do.call(rBind, lapply(1:10, function(i) 
    set.seed(i)
    m[-sample(1:n, n/3), ]
))
m <- t(m) %*% m

一些描述（感谢 Gabor Grothendieck）：

> dim(m)
[1] 100000 100000
> sum(m) / prod(dim(m))
[1] 6.6667e-05
> table(rowSums(m))

    0     1     2     3     4     5     6     7     8     9    10 
    5    28   320  1622  5721 13563 22779 26011 19574  8676  1701 
> table(colSums(m))

    0     1     2     3     4     5     6     7     8     9    10 
    5    28   320  1622  5721 13563 22779 26011 19574  8676  1701 

还有一些错误信息：

> Matrix::solve(m)
Error in LU.dgC(a) : cs_lu(A) failed: near-singular A (or out of memory)
> base::solve(m)
Error in asMethod(object) : 
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105
> MASS::ginv(m)
Error in MASS::ginv(m) : 'X' must be a numeric or complex matrix

赏金的目标是找到一种方法来计算 m 的 Moore-Penrose 广义逆，而 RAM 小于 8Gb（速度和性能并不重要）。

【问题讨论】：

计算出逆后需要做什么？ mathoverflow.net/questions/72616/…

谢谢@BenBolker，我尝试使用大矩阵在scribd.com/doc/226039409#page=14 的第14 页实现公式（22）。对于如何保存该计算，我也将不胜感激，因为使用 Matrix::sparseMatrix 计算 D_1 和 D^a_2 非常快，但我被 (pseduo) 逆运算卡住了。

Scribd 说这是一份私人文件……

@BenBolker 对于这个蹩脚的问题，我感到非常抱歉，请再试一次，我刚刚更新了隐私设置。

由于某种原因，您的示例 m 没有创建方阵。

【参考方案1】：

如果您的 1 很少，那么您的矩阵在任何列和任何行中的 1 都可能不超过一个，在这种情况下，广义逆等于转置：

library(Matrix)
set.seed(123)
n <- 1e5
m <- sparseMatrix(sample(1:n, n/10), sample(1:n, n/10), x = 1, dims = c(n, n))

##############################################################################
# confirm that m has no more than one 1 in each row and column
##############################################################################
table(rowSums(m))  # here 90,000 rows are all zero, 10,000 have a single one
##     0     1 
## 90000 10000 

table(colSums(m))  # here 90,000 cols are all zero, 10,000 have a single one
##     0     1 
## 90000 10000 

##############################################################################
# calculate generalized inverse
##############################################################################
minv <- t(m)

##############################################################################
# check that when multiplied by minv it gives a diagonal matrix of 0's and 1's
##############################################################################
mm <- m %*% minv

table(diag(mm)) # diagonal has 90,000 zeros and 10,000 ones
##     0     1 
## 90000 10000 

diag(mm) <- 0
range(mm) # off diagonals are all zero
## [1] 0 0

修订改进了演示。

【讨论】：

太棒了，非常感谢@G-Grothendieck！不幸的是，m 在某些列/行中有几个1s，尽管它确实很稀疏：sum(m) / prod(dim(m)) 在我的矩阵上返回0.0007337982。我试着想出一个更好的例子。

什么是table(rowSums(m)) 和table(colSums(m)) 和dim(m)？

@G-Grothendieck 请在gist.github.com/daroczig/9bc956ce8e1cb210234a 找到详细信息，直到我想出一个真正的最小可重现示例。

这似乎比我想象的要密集，但你可以像这样减小问题的大小： ix 0), which(colSums(m) > 0)); minv ginv 不起作用，但如果你能找到另一个可能会减少它足以提供帮助。

太棒了，非常感谢@G-Grothendieck！我会将赏金开放几天，但我很确定提供的声望点最终会出现在您的帐户中。再次感谢！