How can I compute the null space/kernel (x: M·x = 0) of a sparse matrix in Python?

Question

I found some examples online showing how to find the null space of a regular matrix in Python, but I couldn't find any examples for a sparse matrix (scipy.sparse.csr_matrix).

By null space I mean x such that M·x = 0, where '·' is matrix multiplication. Does anybody know how to do this?

Furthermore, in my case I know that the null space will consist of a single vector. Can this information be used to improve the efficiency of the method?

Answer 1

This isn't a complete answer yet, but hopefully it will be a starting point towards one. You should be able to compute the null space using a variant on the SVD-based approach shown for dense matrices in this question:

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

If you know that x is a single row vector then in principle you would only need to compute the smallest singular value/vector of M. It ought to be possible to do this using scipy.sparse.linalg.svds, i.e.:

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

Unfortunately, scipy's svds seems to be badly behaved when finding small singular values of singular or near-singular matrices and usually either returns NaNs or throws an ArpackNoConvergence error.

I'm not currently aware of an alternative implementation of truncated SVD with Python bindings that will work on sparse matrices and can selectively find the smallest singular values - perhaps someone else knows of one?

Edit

As a side note, the second approach seems to work reasonably well using MATLAB or Octave's svds function:

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

Answer 2

I have been trying to find a solution to the same problem. Using Scipy's svds function provides unreliable results for small singular values. Therefore i have been using QR decomposition instead. The sparseqr https://github.com/yig/PySPQR provides a wrapper for Matlabs SuiteSparseQR method, and works reasonably well. Using this the null space can be calculated as:

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