python - SymPy rref() 返回奇异矩阵的单位矩阵

Question

import numpy
import sympy

n = 7
k = 3

X = numpy.random.randn(n,k)
Px = X@numpy.linalg.inv(numpy.transpose(X)@X)@numpy.transpose(X) #X(X'X)^(-1)X'

print(sympy.Matrix(Px).rref())

您可以自己验证，Px 是单数。但是，sympy.rref() 返回:

(Matrix([[1, 0, 0, 0, 0, 0, 0],

[0, 1, 0, 0, 0, 0, 0],

[0, 0, 1, 0, 0, 0, 0],

[0, 0, 0, 1, 0, 0, 0],

[0, 0, 0, 0, 1, 0, 0],

[0, 0, 0, 0, 0, 1, 0],

[0, 0, 0, 0, 0, 0, 1]]), (0, 1, 2, 3, 4, 5, 6))

为什么它不返回真正的 rref？我在某处读到可以传递 simple=True，但是它没有任何区别。

Answer 1

In [49]: Px                                                                     
Out[49]: 
array([[ 0.5418898 ,  0.44245552,  0.04973693, -0.06834885, -0.19086119,
        -0.07003176,  0.06325021],...
       [ 0.06325021, -0.11080081,  0.21656224, -0.07445145, -0.28634725,
         0.06648907,  0.19199866]])
In [50]: np.linalg.det(Px)                                                      
Out[50]: 2.141647537907433e-67
In [51]: np.linalg.inv(Px)                                                      
Out[51]: 
array([[-7.18788695e+15,  4.95655702e+15,  7.52738018e+15,
        -4.40875311e+15, -1.64015565e+16,  2.63785320e+15,
        -3.03465003e+16],
       [ 1.59176426e+16, ....
       [ 3.31636798e+16, -3.39094560e+16, -3.60287970e+16,
        -1.27160460e+16,  2.14338015e+16,  3.32345350e+15,
         3.60287970e+16]])

您的 Px 接近单数，但又不完全如此。对比一下

In [52]: M = np.arange(9).reshape(3,3)                                          
In [53]: np.linalg.det(M)                                                       
Out[53]: 0.0
In [55]: np.linalg.inv(M)                                                       
LinAlgError: Singular matrix

In [56]: sympy.Matrix(M).rref()                                                 
Out[56]: 
(Matrix([
 [1, 0, -1],
 [0, 1,  2],
 [0, 0,  0]]), (0, 1))

从数字上讲，您的 Px 不是单数，只是接近:

In [57]: sympy.Matrix(Px).rref()                                                
Out[57]: 
(Matrix([
 [1, 0, 0, 0, 0, 0, 0],
 [0, 1, 0, 0, 0, 0, 0],
 [0, 0, 1, 0, 0, 0, 0],
 [0, 0, 0, 1, 0, 0, 0],
 [0, 0, 0, 0, 1, 0, 0],
 [0, 0, 0, 0, 0, 1, 0],
 [0, 0, 0, 0, 0, 0, 1]]), (0, 1, 2, 3, 4, 5, 6))

但是使用自定义的 iszerofunc:

In [58]: sympy.Matrix(Px).rref(iszerofunc=lambda x: abs(x)<1e-16)               
Out[58]: 
(Matrix([
 [1, 0, 0,  0.647383887198708, -1.91409951634531, -1.43377991000974,  0.578981680134581],
 [0, 1, 0, -0.839184067893959,  1.88998490600173,  1.43367640627271, -0.611620902311026],
 [0, 0, 1, -0.962221703397948, 0.203783478612254,  1.45929622452135,  0.404548167005728],
 [0, 0, 0,                  0,                 0,                 0,                  0],
 [0, 0, 0,                  0,                 0,                 0,                  0],
 [0, 0, 0,                  0,                 0,                 0,                  0],
 [0, 0, 0,                  0,                 0,                 0,                  0]]),
 (0, 1, 2))