如何在 Python 中计算稀疏矩阵的零空间/内核 (x: M·x = 0)？答案

【问题标题】：How can I compute the null space/kernel (x: M·x = 0) of a sparse matrix in Python?如何在 Python 中计算稀疏矩阵的零空间/内核 (x: M·x = 0)？
【发布时间】：2016-01-29 08:29:12
【问题描述】：

我在网上找到了一些示例，展示了如何在 Python 中找到正则矩阵的零空间，但我找不到任何稀疏矩阵的示例 (scipy.sparse.csr_matrix)。

空空间是指 x 使得 M·x = 0，其中 ' ·' 是矩阵乘法。有人知道怎么做这个吗？

此外，就我而言，我知道零空间将由单个向量组成。这些信息可以用来提高方法的效率吗？

【问题讨论】：

我不太清楚你的意思。简单地说，M * x == 0 对于全零向量来说是正确的，但我认为这不是你想要的。
我想你一定是指矩阵的row null space，即当你说 "regular row/column multiplication" 你指的是矩阵向量而不是比元素乘法。
scipy.sparse.linalg.spsolve 有点这样做，但我猜它不如直接查找空空间有效。根据个人经验（b 非零），它非常慢。
您好，感谢您的快速回复！是的，我指的是矩阵向量乘法，我正在寻找一个非零元素 x 使得 M * x = 0。这在行列式 = 0 的矩阵中是可能的，或者如果矩阵的形状不是正方形。谢谢 simonzack，实际上这正是我目前正在做的事情，但我想知道是否有更有效的方法来找到空空间，因为 spsolve 产生了更通用的答案。
你能给我们举一个密集数组的小例子吗？在考虑稀疏等价物时，这将作为参考点。

标签： python numpy scipy linear-algebra sparse-matrix

【解决方案1】：

这还不是一个完整的答案，但希望这将是一个起点。您应该能够使用针对密集矩阵in this question 显示的基于 SVD 的方法的变体计算零空间：

import numpy as np
from scipy import sparse
import scipy.sparse.linalg


def rand_rank_k(n, k, **kwargs):
    "generate a random (n, n) sparse matrix of rank <= k"
    a = sparse.rand(n, k, **kwargs)
    b = sparse.rand(k, n, **kwargs)
    return a.dot(b)

# I couldn't think of a simple way to generate a random sparse matrix with known
# rank, so I'm currently using a dense matrix for proof of concept
n = 100
M = rand_rank_k(n, n - 1, density=1)

# # this seems like it ought to work, but it doesn't
# u, s, vh = sparse.linalg.svds(M, k=1, which='SM')

# this works OK, but obviously converting your matrix to dense and computing all
# of the singular values/vectors is probably not feasible for large sparse matrices
u, s, vh = np.linalg.svd(M.todense(), full_matrices=False)

tol = np.finfo(M.dtype).eps * M.nnz
null_space = vh.compress(s <= tol, axis=0).conj().T

print(null_space.shape)
# (100, 1)
print(np.allclose(M.dot(null_space), 0))
# True

如果你知道 x 是一个单行向量，那么原则上你只需要计算 M 的最小奇异值/向量。使用scipy.sparse.linalg.svds 应该可以做到这一点，即：

u, s, vh = sparse.linalg.svds(M, k=1, which='SM')
null_space = vh.conj().ravel()

不幸的是，scipy 的 svds seems to be badly behaved 在查找奇异或接近奇异矩阵的小奇异值时，通常会返回 NaN 或抛出 ArpackNoConvergence 错误。

我目前不知道使用 Python 绑定的截断 SVD 的替代实现，它可以在稀疏矩阵上工作并且可以选择性地找到最小的奇异值 - 也许其他人知道？

编辑

附带说明一下，使用 MATLAB 或 Octave 的 svds 函数，第二种方法似乎工作得相当好：

>> M = rand(100, 99) * rand(99, 100);
% svds converges much more reliably if you set sigma to something small but nonzero
>> [U, S, V] = svds(M, 1, 1E-9);
>> max(abs(M * V))
ans =    1.5293e-10

【讨论】：

【解决方案2】：

我一直在尝试找到解决同一问题的方法。使用 Scipy 的 svds 函数会为小的奇异值提供不可靠的结果。因此，我一直在使用 QR 分解。 sparseqr https://github.com/yig/PySPQR 为 Matlabs SuiteSparseQR 方法提供了一个包装器，并且运行良好。使用这个零空间可以计算为：

    from sparseqr import qr
    Q, _, _,r = qr( M.transpose() )
    N = Q.tocsr()[:,r:]

【讨论】：