Python (NumPy, SciPy), finding the null space of a matrix

Question

I'm trying to find the null space (solution space of Ax=0) of a given matrix. I've found two examples, but I can't seem to get either to work. Moreover, I can't understand what they're doing to get there, so I can't debug. I'm hoping someone might be able to walk me through this.

The documentation pages (numpy.linalg.svd, and numpy.compress) are opaque to me. I learned to do this by creating the matrix C = [A|0], finding the reduced row echelon form and solving for variables by row. I can't seem to follow how it's being done in these examples.

Thanks for any and all help!

Here is my sample matrix, which is the same as the wikipedia example:

A = matrix([
    [2,3,5],
    [-4,2,3]
    ])  

Method (found here, and here):

import scipy
from scipy import linalg, matrix
def null(A, eps=1e-15):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

When I try it, I get back an empty matrix:

Python 2.6.6 (r266:84292, Sep 15 2010, 16:22:56) 
[GCC 4.4.5] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import scipy
>>> from scipy import linalg, matrix
>>> def null(A, eps=1e-15):
...    u, s, vh = scipy.linalg.svd(A)
...    null_mask = (s <= eps)
...    null_space = scipy.compress(null_mask, vh, axis=0)
...    return scipy.transpose(null_space)
... 
>>> A = matrix([
...     [2,3,5],
...     [-4,2,3]
...     ])  
>>> 
>>> null(A)
array([], shape=(3, 0), dtype=float64)
>>> 

Solution 1

It appears to be working okay for me:

A = matrix([[2,3,5],[-4,2,3],[0,0,0]])
A * null(A)
>>> [[  4.02455846e-16]
>>>  [  1.94289029e-16]
>>>  [  0.00000000e+00]]

Solution 2

Sympy makes this straightforward.

>>> from sympy import Matrix
>>> A = [[2, 3, 5], [-4, 2, 3], [0, 0, 0]]
>>> A = Matrix(A)
>>> A * A.nullspace()[0]
Matrix([
[0],
[0],
[0]])
>>> A.nullspace()
[Matrix([
[-1/16],
[-13/8],
[    1]])]

Solution 3

As of last year (2017), scipy now has a built-in null_space method in the scipy.linalg module (docs).

The implementation follows the canonical SVD decomposition and is pretty small if you have an older version of scipy and need to implement it yourself (see below). However, if you're up-to-date, it's there for you.

def null_space(A, rcond=None):
    u, s, vh = svd(A, full_matrices=True)
    M, N = u.shape[0], vh.shape[1]
    if rcond is None:
        rcond = numpy.finfo(s.dtype).eps * max(M, N)
    tol = numpy.amax(s) * rcond
    num = numpy.sum(s > tol, dtype=int)
    Q = vh[num:,:].T.conj()
    return Q

Solution 4

You get the SVD decomposition of the matrix A. s is a vector of eigenvalues. You are interested in almost zero eigenvalues (see $A*x=\lambda*x$ where $\abs(\lambda)<\epsilon$), which is given by the vector of logical values null_mask.

Then, you extract from the list vh the eigenvectors corresponding to the almost zero eigenvalues, which is exactly what you are looking for: a way to span the null space. Basically, you extract the rows and then transpose the results so that you get a matrix with eigenvectors as columns.

Solution 5

A faster but less numerically stable method is to use a rank-revealing QR decomposition, such as scipy.linalg.qr with pivoting=True:

import numpy as np
from scipy.linalg import qr

def qr_null(A, tol=None):
    Q, R, P = qr(A.T, mode='full', pivoting=True)
    tol = np.finfo(R.dtype).eps if tol is None else tol
    rnk = min(A.shape) - np.abs(np.diag(R))[::-1].searchsorted(tol)
    return Q[:, rnk:].conj()

For example:

A = np.array([[ 2, 3, 5],
              [-4, 2, 3],
              [ 0, 0, 0]])
Z = qr_null(A)

print(A.dot(Z))
#[[  4.44089210e-16]
# [  6.66133815e-16]
# [  0.00000000e+00]]

View more solutions

Comments

[Valeria]

I'm sure I'm missing something, but Wikipedia suggests that the values should be [ [-.0625], [-1.625], [1] ]?

[Valeria]

Moreover, it's returning an empty matrix for me []. What could be wrong?

[Joe Kington]

@Nona Urbiz - It's returning an empty matrix because you're not putting in a row of zeros, as Bashwork (and wikipedia) does above. Also, the null space values returned ([-0.33, -0.85, 0.52]) are normalized so that the magnitude of the vector is 1. The wikipedia example is not normalized. If you just take n = null(A) and have a look at n / n.max(), you'll get [-.0625, -1.625, 1].

[Valeria]

Thank you very much for taking the time to reply and help me. Your answer was very helpful to me, but I accepted Bashworks answer as ultimately, it brought me to the solution. The only reason I am able to understand the solution though, is your response.

[Wok]

No worry, I thought your problem was something else.

[Coder]

@Bashwork - How would I know to programmatically add a row of zeroes? Does the matrix have to be square?

[Thomas Wentworth]

I have been using this code in my work and I noticed a problem. An eps value of 1e-15 seems to be too small. Notably, consider the matrix A = np.ones(13,2). This code will report that this matrix has a rank 0 null space. This is due to the scipy.linalg.svd function reporting that the second singular value is above 1e-15. I don't know much about the algorithms behind this function, however I suggest using eps=1e-12 (and perhaps lower for very large matrices) unless someone with more knowledge can chime in. (In infinite precision the second singular value should be 0).

[cxxl]

Exactly the solution I was looking for!

[Adriaan]

I believe this approach only works for integers but not for floats :(.

[Grant Curell]

I tried this out and got different answers on Matrix([[2, 2, 1, 0], [2, -2, -1, 1]]). Using sympy I correctly got: ``` [Matrix([ [ 0], [-1/2], [ 1], [ 0]]), Matrix([ [-1/4], [ 1/4], [ 0], [ 1]])] ``` But with the above method I got: ``` array([[0.00000000e+00, 4.16333634e-17], [1.73472348e-16, 0.00000000e+00]]) ``` Am I missing something? (I'm decent at linear algebra but I have no idea what this method is - I'm just trying it)

[Grant Curell]

@Adriaan I just tested with floats and it worked for me.