使用 numpy 计算零空间的有理基础

Question

我正在尝试计算矩阵零空间的有理基础。有很多关于如何使用 Python/numpy 计算零空间的帖子，但他们计算的是正交基而不是有理基。这是在 MATLAB 中如何完成的：

ns = null(A,'r')

看源码的时候看到是这样计算的：

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

这里rref 是简化的行梯形。因此，通过查看此源代码，我尝试使用以下代码重新创建它：

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

有两件事我不知道。首先，我不知道如何将比率部分添加到rref 函数中。其次，我不确定我的索引是否正确，因为 MATLAB 的索引从 1 开始，并且当您选择切片时索引包括最后一个元素（即 1:5 包括 1 和 5）。

Answer 1

SymPy 开箱即用，尽管（在 Python 中是象征性的）不如 NumPy 或 Scipy 快。浮点输入示例：

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

打印两个列向量的列表，表示为一列矩阵。

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

nsimplify 的使用对于将浮点数转换为它们要表示的有理数是必要的。如果矩阵创建为整数/有理项矩阵，则没有必要。

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

这里使用S(3)/2 代替`3/2 以强制创建SymPy 对象而不是浮点计算。