Calculating rational basis for the nullspace using numpy

Question

I am trying to calculate the rational basis for null space of a matrix. There is quite a few posts about how nullspace is calculated using Python/numpy but they calculate it for orthonormal basis and not for the rational basis. Here is how this is done in MATLAB:

ns = null(A,'r')

When I look at the source code, I saw that it is calculated like this:

function Z = null(A,how)
[m,n] = size(A)
%...
[R,pivcol] = rref(A);
r = length(pivcol);
nopiv = 1:n;
nopiv(pivcol) = [];
Z = zeros(n,n-r,class(A));
if n > r
  Z(nopiv,:) = eye(n-r,n-r,class(A));
  if r > 0
     Z(pivcol,:) = -R(1:r,nopiv);
  end
end
%...

function [A,jb] = rref(A,tol)
%...
[m,n] = size(A);

[num, den] = rat(A);
rats = isequal(A,num./den);

if (nargin < 2), tol = max(m,n)*eps(class(A))*norm(A,'inf'); end

i = 1;
j = 1;
jb = [];
while (i <= m) && (j <= n)

   [p,k] = max(abs(A(i:m,j))); k = k+i-1;
   if (p <= tol)

      A(i:m,j) = zeros(m-i+1,1);
      j = j + 1;
   else

      jb = [jb j];

      A([i k],j:n) = A([k i],j:n);

      A(i,j:n) = A(i,j:n)/A(i,j);

      for k = [1:i-1 i+1:m]
         A(k,j:n) = A(k,j:n) - A(k,j)*A(i,j:n);
      end
      i = i + 1;
      j = j + 1;
   end
end

if rats
    [num,den] = rat(A);
    A=num./den;
end

Here rref is the reduced row echelon form. Thus by looking at this source code I tried to recreate it with following code:

def fract(x):
    return Fraction(x)

def dnm(x):
    return x.denominator

def nmr(x):
    return x.numerator

fractionize = np.vectorize(fract)
denom = np.vectorize(dnm)
numer = np.vectorize(nmr)

def rref(A,tol=1e-12):

    m,n = A.shape
    Ar = A.copy()
    i,j = 0,0
    jb = []
    while i < m and j < n:
        p = np.max(np.abs(Ar[i:m,j]))
        k = np.where(np.abs(Ar[i:m,j]) == p)[0][0]
        k = k + i - 1
        if (p <= tol):
            Ar[i:m,j] = np.zeros((m-i,))
            j += 1
        else:
            jb.append(j)
            Ar[(i,k),j:n] = Ar[(k,i),j:n]
            Ar[i,j:n] = Ar[i,j:n]/Ar[i,j]
            for k in np.hstack((np.arange(0,i),np.arange(i+1,m))):
                Ar[k,j:n] = Ar[k,j:n] - Ar[k,j]*A[i,j:n]
            i += 1
            j += 1
    print(len(jb))
    return Ar,jb

def null(A,tol=1e-5):

    m,n = A.shape
    R,pivcol = rref(A,tol=tol)
    print(pivcol)
    r = len(pivcol)
    nopiv = np.ones(n).astype(bool)
    nopiv[pivcol] = np.zeros(r).astype(bool)
    Z = np.zeros((n,n-r))
    if n > r:
        Z[nopiv,:] = np.eye(n-r,n-r)
        if r > 0:
            Z[pivcol,:] = -R[:r,nopiv]
    return Z

There are two things that I don't know. First, I do not know how to add the ratios part into rref function. Second, I am not sure if my indexes are correct since MATLAB's indices are start from 1 and indexing includes the last element when you choose for a slice (i.e. 1:5 includes both 1 and 5).

Answer 1

SymPy does that out of the box, although (being symbolic, and in Python) not as fast as NumPy or Scipy would. An example with floating point input:

from sympy import Matrix, S, nsimplify
M = Matrix([[2.75, -1.2, 0, 3.2], [8.29, -4.8, 7, 0.01]])
print(nsimplify(M, rational=True).nullspace())

Prints a list of two column vectors, represented as one-column matrices.

[Matrix([
[  700/271],
[9625/1626],
[        1],
[        0]]), Matrix([
[  -1279/271],
[-17667/2168],
[          0],
[          1]])]

The use of nsimplify was necessary to convert floats to the rationals that they were meant to represent. If the matrix is created as a matrix of integer/rational entries, that would not be necessary.

M = Matrix([[1, 2, 3, 5, 9], [9, -3, 0, 2, 4], [S(3)/2, 0, -1, 2, 0]])
print(M.nullspace())

[Matrix([
[ -74/69],
[-176/69],
[   9/23],
[      1],
[      0]]), Matrix([
[ -70/69],
[-118/69],
[ -35/23],
[      0],
[      1]])]

Here, S(3)/2 is used instead of `3/2 in order to force SymPy object creation instead of floating point evaluation.