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

Question

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

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

此外,在我的情况下,我知道零空间将由单个向量组成.这些信息可以用来提高方法的效率吗？

Answer 1

这还不是一个完整的答案，但希望它将成为一个起点。您应该能够使用本问题中针对稠密矩阵所示的基于 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()

不幸的是，当寻找奇异或接近奇异矩阵的小奇异值时，scipysvds 似乎表现得很糟糕，并且通常返回 NaN 或抛出错误ArpackNoConvergence。

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

编辑

附带说明一下，第二种方法似乎使用 MATLAB 或 Octavesvds函数效果相当好：

>> 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