CODEWITHSUNDEEP
numpymatrixscipylinear-algebrasvd
Nona Urbiz
Nona Urbiz

Reputation: 5013

Finding the null space of a matrix

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. Can someone 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)
>>>

Upvotes: 39

Views: 67062

Answers (7)

Glen
Glen

Reputation: 1

import numpy as np
from scipy.linalg import null_space
A = np.random.randn(5,8)
n = null_space(A)

Upvotes: 0

sputnick
sputnick

Reputation: 323

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

Upvotes: 16

Idr
Idr

Reputation: 6250

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

Upvotes: 27

ali_m
ali_m

Reputation: 74232

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

Upvotes: 4

Thomas Wentworth
Thomas Wentworth

Reputation: 41

Your method is almost correct. The issue is that the shape of s returned by the function scipy.linalg.svd is (K,) where K=min(M,N). Thus, in your example, s only has two entries (the singular values of the first two singular vectors). The following correction to your null function should allow it to work for any sized matrix.

import scipy
import numpy as np
from scipy import linalg, matrix
def null(A, eps=1e-12):
...    u, s, vh = scipy.linalg.svd(A)
...    padding = max(0,np.shape(A)[1]-np.shape(s)[0])
...    null_mask = np.concatenate(((s <= eps), np.ones((padding,),dtype=bool)),axis=0)
...    null_space = scipy.compress(null_mask, vh, axis=0)
...    return scipy.transpose(null_space)
A = matrix([[2,3,5],[-4,2,3]])
print A*null(A)
>>>[[  4.44089210e-16]
>>> [  6.66133815e-16]]
A = matrix([[1,0,1,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]])
print null(A)
>>>[[ 0.         -0.70710678]
>>> [ 0.          0.        ]
>>> [ 0.          0.70710678]
>>> [ 1.          0.        ]]
print A*null(A)
>>>[[ 0.  0.]
>>> [ 0.  0.]
>>> [ 0.  0.]
>>> [ 0.  0.]]

Upvotes: 3

Bashwork
Bashwork

Reputation: 1619

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

Upvotes: 7

Wok
Wok

Reputation: 5333

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.

Upvotes: 13

Related Questions