如何计算大型稀疏矩阵的伪逆？解决CSR矩阵求逆函数缺失问题

针对你遇到的50000×12000大型稀疏矩阵求伪逆的问题，以下是几种可行的方案：

1. 基于稀疏SVD构造伪逆 #

稀疏矩阵的伪逆可通过奇异值分解（SVD）推导：若矩阵$A$的SVD为$A = U\Sigma V^{T$，则伪逆$A}+ = V\Sigma^+ U^{T$，其中$\Sigma}+$是$\Sigma$的逆矩阵（仅保留非零奇异值的倒数，其余为0）。

对于大型稀疏矩阵，无需计算全量SVD，使用scipy.sparse.linalg.svds计算部分重要的奇异值和奇异向量即可，大幅节省内存：

import scipy.sparse.linalg as spla
import scipy.sparse as sp
import numpy as np

# 假设A是已转换好的csr_matrix
# k取前k个最大的奇异值，可根据方差解释率或矩阵估计秩调整，比如设为1000
u, s, vt = spla.svds(A, k=1000)

# 构造奇异值的伪逆对角矩阵
s_inv = np.diag(1 / s)
# 计算伪逆
A_pinv = vt.T @ s_inv @ u.T
# 若需要转回稀疏格式
A_pinv_sparse = sp.csr_matrix(A_pinv)

2. 利用超定系统的伪逆公式 #

由于你的矩阵是超定的（行数50000 > 列数12000），伪逆可表示为$A^+ = (A^T A)^{-1} A^T$。$AT A$是12000×12000的矩阵，尺寸远小于原矩阵，计算压力小：

import numpy as np
import scipy.sparse as sp

# 计算A^T A
ata = A.T @ A
# 若A列满秩，直接求逆；否则用伪逆
if np.linalg.matrix_rank(ata.toarray()) == ata.shape[0]:
    ata_inv = np.linalg.inv(ata.toarray())
else:
    ata_inv = np.linalg.pinv(ata.toarray())

# 计算伪逆
A_pinv = ata_inv @ A.T.toarray()
# 转稀疏格式（可选）
A_pinv_sparse = sp.csr_matrix(A_pinv)

3. 迭代法避免显式计算伪逆 #

如果你的最终需求是计算$A^+ b$（伪逆作用于某个向量），无需显式构造伪逆矩阵，直接用迭代最小二乘求解器即可，内存效率最高：

import scipy.sparse.linalg as spla

# 求解Ax = b的最小二乘解，等价于x = A^+ b
x, info = spla.lsqr(A, b)

方案选择建议 #

• 若需要显式伪逆矩阵：优先选择稀疏SVD方案（适合低秩矩阵），或超定公式方案（适合列满秩矩阵）。
• 若仅需伪逆的作用结果：优先用迭代法，避免存储庞大的伪逆矩阵。

内容的提问来源于stack exchange，提问作者Hold My Stack