使用numpy计算无效空间的有理基础 [英] Calculating rational basis for the nullspace using numpy
问题描述
我正在尝试计算矩阵零空间的有理基础.关于如何使用Python/numpy计算零空间的文章很多，但是他们是基于正交的基础而不是有理基础的.这是在MATLAB中完成的方法:
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
此处rref是精简行梯形形式.因此，通过查看此源代码，我尝试使用以下代码重新创建它:
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
我不知道有两件事.首先，我不知道如何将比率部分添加到rref函数中.其次，我不确定我的索引是否正确，因为MATLAB的索引从1开始，并且在选择切片时索引包括最后一个元素(即1:5包括1和5).
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).
推荐答案
SymPy 开箱即用，尽管(作为符号，在Python中)不如NumPy或Scipy快.带有浮点输入的示例:
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]])]
必须使用nsimplify才能将浮点数转换为它们要表示的基本原理.如果将矩阵创建为整数/有理数条目的矩阵，则没有必要.
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]])]
此处，使用S(3)/2代替`3/2来强制创建SymPy对象而不是浮点运算.
Here, S(3)/2 is used instead of `3/2 in order to force SymPy object creation instead of floating point evaluation.
这篇关于使用numpy计算无效空间的有理基础的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
