为何Numpy中matmul与linalg.solve求解线性方程组结果不一致

两种求解方法结果差异的根本原因 #

你所用的系数矩阵A是奇异（不可逆）矩阵，不存在严格的数学逆，两种方法的数值稳定性差异直接导致了结果不同：

• 矩阵A的三行线性相关：第二行减第一行等于[3,3,3]，第三行减第二行也等于[3,3,3]，其数学意义上的行列式为0，不满足矩阵可逆的充要条件。此时调用np.linalg.inv(A)得到的结果本身就是浮点误差下的无效输出，用这个错误的逆矩阵和b做乘法np.matmul(np.linalg.inv(A), b)，自然得不到合理的解。
• np.linalg.solve不会显式计算矩阵逆，它会通过LU分解等数值稳定性更高的算法直接求解线性方程组，对于这类不满秩的方阵，会返回最小二乘意义下残差最小的近似解，你可以验证将solve返回的结果代入A @ x，得到的输出和b的误差在浮点精度范围内，是符合方程组近似解要求的。
• 补充说明：对于不可逆的系数矩阵，求解线性方程组的标准做法是使用最小二乘函数np.linalg.lstsq，而不是显式求逆或者直接调用solve，可以进一步提升结果的可靠性。

内容的提问来源于stack exchange，提问作者Anne Maier