SymPy的秩与NumPy矩阵秩不同

4

给定一个 SymPy 矩阵 M

M = Matrix([
 [0.000111334436666596, 0.00114870370895408, -0.000328330524152990, 5.61388353859808e-6, -0.000464532588930332, -0.000969955779635878, 1.70579589853818e-5, -5.77891177019884e-6, -0.000186812539472235, -2.37115911398055e-5],
 [-0.00105346453420510, 0.000165063406707273, -0.00184449574409890, 0.000658080565333929, 0.00197652092300241, 0.000516180213512589, 9.53823860082390e-5, 0.000189858427211978, -3.80494288487685e-5, 0.000188984043643408],
 [-0.00102465075104153, -0.000402915220398109, 0.00123785300884241, -0.00125808154543978, 0.000126618511490838, 0.00185985865307693, 0.000123626008509804, 0.000211557638637554, 0.000407232404255796, 1.89851719447102e-5],
 [0.230813497584639, -0.209574389008468, 0.742275067362657, -0.202368828927654, -0.236683258718819, 0.183258819107153, 0.180335891933511, -0.530606389541138, -0.379368598768419, 0.334800403899511],
 [-0.00102465075104153, -0.000402915220398109, 0.00123785300884241, -0.00125808154543978, 0.000126618511490838, 0.00185985865307693, 0.000123626008509804, 0.000211557638637554, 0.000407232404255796, 1.89851719447102e-5],
 [0.00105346453420510, -0.000165063406707273, 0.00184449574409890, -0.000658080565333929, -0.00197652092300241, -0.000516180213512589, -9.53823860082390e-5, -0.000189858427211978, 3.80494288487685e-5, -0.000188984043643408],
 [0.945967255845168, -0.0468645728473480, 0.165423896937049, -0.893045423193559, -0.519428986944650, -0.0463256408085840, -0.0257001217930424, 0.0757328764368606, 0.0541336731317414, -0.0477734271777646],
 [-0.0273371493900004, -0.954100482348723, -0.0879282784854250, 0.100704543595514, -0.243312734473589, -0.0217088779350294, 0.900584332231093, 0.616061129532614, 0.0651163853434486, -0.0396603397583054],
 [0.0967584768347089, -0.0877680087304911, -0.667679934757176, -0.0848411039101494, -0.0224646387789634, -0.194501966574153, 0.0755161040544943, 0.699388977592066, 0.394125039254254, -0.342798611994521],
 [-0.000222668873333193, -0.00229740741790816, 0.000656661048305981, -1.12277670771962e-5, 0.000929065177860663, 0.00193991155927176, -3.41159179707635e-5, 1.15578235403977e-5, 0.000373625078944470, 4.74231822796110e-5]
 ])

我已经用 SymPy 计算了矩阵的 rank()rref()。其中,它的秩是 7rref() 的结果如下:
Matrix([
[1, 0, 0, 0, 0, 0, 0, -5.14556976678473, -3.72094268951566,  3.48581267477014],
[0, 1, 0, 0, 0, 0, 0, -5.52930150663022, -4.02230308325653,  3.79193678096199],
[0, 0, 1, 0, 0, 0, 0,  2.44893308665325,  1.83777402439421, -1.87489784909824],
[0, 0, 0, 1, 0, 0, 0, -7.33732284392352, -5.25036238623229,  4.97256759287563],
[0, 0, 0, 0, 1, 0, 0,  5.48049237370489,  3.90091366576548, -3.83642187384021],
[0, 0, 0, 0, 0, 1, 0, -10.6826798792866, -7.56560803870182,  7.45974067056387],
[0, 0, 0, 0, 0, 0, 1, -3.04726210012149, -2.66388837034592,  2.48327234504403],
[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,                 0]])

奇怪的事情是,如果我使用NumPyMATLAB计算排名,我得到的值为6,如果我用MATLAB计算rref,我得到了预期的结果——最后4行都是零(而不仅仅是最后3行)。

有没有人知道这种差异来自哪里,以及为什么我无法使用SymPy获得正确的结果?我知道排名6是正确的,因为它是一个存在某些线性依赖关系的方程组。

1个回答

3

从矩阵的特征值来看,排名确实为6:

array([ 1.14550481e+00+0.00000000e+00j, -1.82137718e-01+6.83443168e-01j,
   -1.82137718e-01-6.83443168e-01j,  2.76223053e-03+0.00000000e+00j,
   -3.51138883e-04+8.61508469e-04j, -3.51138883e-04-8.61508469e-04j,
    5.21160131e-17+0.00000000e+00j, -2.65160469e-16+0.00000000e+00j,
   -2.67753616e-18+9.70937977e-18j, -2.67753616e-18-9.70937977e-18j])

使用我拥有的sympy版本,与numpy返回的6级相比,我甚至得到了8级的排名。

但实际上,由于矩阵的大小(可能与SymPy无法计算此矩阵的特征值有关),Sympy无法解决该矩阵的特征值。

因此,其中一个Sympy试图通过符号方程求解并找到排名(基于不完美的浮点数),而另一个numpy使用近似值(lapack IIRC)来查找特征值。通过具有适当阈值,numpy找到了正确的排名,但在不同的阈值下可能会得出不同的结果。 Sympy尝试基于完美的6级系统的近似系统来找到排名,并发现它是7或8级。由于浮点差异,这并不奇怪(例如,Sympy移动到整数以尝试找到特征值,而不是停留在浮点领域中)。


我现在明白问题来自哪里了。您是否知道可能的解决方案? 我尝试过从互联网上调整/修改现有的Python rref函数,但到目前为止没有找到真正的解决方案。 - carobnodrvo
获取排名的目标是什么?不幸的是,在浮点情况下,行/列并不真正共线,这使得符号解决方案难以使用。 - Matthieu Brucher
嗯,实际上我并不是在计算排名,而是要减少呈现矩阵(以及其他类似的矩阵)的行数。 - carobnodrvo

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接