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

Question

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

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

我尝试过的:

• solve 产生 LU.dgC(a) 中的错误:cs_lu(A) 失败:接近奇异 A(或内存不足) 错误
• MASS::ginv 与 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 错误“内存不足”。/Core/cholmod_memory.c，第 147 行错误消息

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

<小时/>

快速可重现的示例:

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

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

Answer 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

修订改进了演示。