9

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?

ali_m
  • 71,714
  • 23
  • 223
  • 298
Paolo Erdman
  • 331
  • 1
  • 6
  • I'm not quite sure what you mean. Trivially, `M * x == 0` would be true for a vector of all zeros, but I assume that's not what you want. – ali_m Oct 29 '15 at 09:40
  • I suppose you must be referring to the [row null space](https://en.wikipedia.org/wiki/Kernel_(linear_algebra)#The_row_space_of_a_matrix) of the matrix, i.e. when you say *"regular row/column multiplication"* you are referring to to matrix-vector rather than elementwise multiplication. – ali_m Oct 29 '15 at 10:24
  • `scipy.sparse.linalg.spsolve` kind of does this, but I'm guessing it's not as efficient as finding the null space directly. From personal experience (b is non-zero) it is quite slow. – simonzack Oct 29 '15 at 10:43
  • Hi, thanks for the quick response! Yes, i am referring the the matrix-vector multiplication, and i am looking for a non-zero element x such that M * x = 0. This is possible in matrices with determinant = 0, or if the shape of the matrix is not square. Thank you simonzack, actually that is exactly what i am currently doing, but i was wondering if there was a more efficient way of finding the null space, since spsolve produces a much more general answear. – Paolo Erdman Oct 29 '15 at 12:10
  • Could you give us a small example with a dense array? That would give as reference point when thinking about a sparse equivalent. – hpaulj Oct 29 '15 at 19:39

2 Answers2

2

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
Community
  • 1
  • 1
ali_m
  • 71,714
  • 23
  • 223
  • 298
1

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:]
Tim Tørnes Pedersen
  • 11
  • 3