Numpy的矩阵解决方案:没有解决方案？

Question

在Numpy中编写一个简单的脚本来解决一个由3个方程组成的系统:

from numpy import *
a = matrix('1 4 1; 4 13 7; 7 22 13')
b = matrix('0;0;1')
print linalg.solve(a,b)

但当我通过命令提示符运行它时,我得到了:

[[  3.46430741e+15]
 [ -6.92861481e+14]
 [ -6.92861481e+14]]

虽然沃尔夫勒姆阿尔法说没有解决方案.

有谁知道为什么Numpy/WRA答案之间似乎存在这种差异？

Answer 1

如果你拿笔和纸(或使用numpy.linalg.matrix_rank)并计算系数和增强矩阵的等级,你会看到,根据Rouché-Capelli定理,没有解决方案.所以Wolfram是对的.

NumPy使用LU分解搜索方程的数值解.没有详细说明LU分解涉及划分,FP算术中的划分可能导致重大错误.

如果你检查:

a = np.matrix('1 4 1; 4 13 7; 7 22 13')
b = np.matrix('0;0;1')

c = np.linalg.solve(a,b)
np.linalg.norm(a * c - b)
## 2.0039024427351748

你会发现NumPy解决方案远非正确.