Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

a simple, matlab-like way of finding the null space of a small matrix in numpy (and number formatting) [duplicate]

There must be a simple way to get a null space of a small (say 3x3) matrix in python's numpy or scipy.

MATLAB can be good about this. Let's say:

A = [1 2 3; 
     2 3 4; 
     2 4 6]

rank(A) % rank is 2 
null(A, 'r') % ask matlab to be ('r') reasonable about 
             % its choice of a vector in A's nullspace

and the output of the last command is:

 1 
-2 
 1

It appears - and is this true? - that things arn't quite as simple in numpy:

import numpy as np
A = array(([1, 2, 3], [2, 3, 4], [2, 4, 6])) 
np.linalg.matrix_rank(A) # ok, getting the rank of a matrix is this esay, even if
                         # it takes more keystrokes, but how about its null space

From what I've searched so far, it appears that one needs to call the svd decomposition function first to get the nullspace.

There must be a simpler way to do this in python.

Also, in matlab one could say:

format rat

to avoid stairing at long decimal fractions in output matrices. (e.g. when the format is set to 'rational' an entry in an output matrix would look like 1/2 instead of the uglier looking 0.5000000)

Does python have a similar feature, or is anyone using python doomed to look at these decimals forever?

Thanks in advance.

d.

like image 507
dave Avatar asked Nov 06 '13 19:11

dave


1 Answers

As you can see from the matlab code for null.m, they also call svd to get the null space. In fact they also call svd to get the rank (I didn't copy that here, but feel free to read the code it's about 4 lines).

function Z = null(A,how)
   % clip
   [~,S,V] = svd(A,0);
   if m > 1, s = diag(S);
      elseif m == 1, s = S(1);
      else s = 0;
   end
   tol = max(m,n) * max(s) * eps(class(A));
   r = sum(s > tol);
   Z = V(:,r+1:n);
end

Here is a python version, this function will compute the rank and null space of a square matrix. It returns the rank and a (N, N - R) array where N is the size of the matrix and R is the rank.

import numpy as np

def null(a, rtol=1e-5):
    u, s, v = np.linalg.svd(a)
    rank = (s > rtol*s[0]).sum()
    return rank, v[rank:].T.copy()
like image 57
Bi Rico Avatar answered Oct 30 '22 11:10

Bi Rico