为什么 Mathematica 和 Python 的答案在处理奇异矩阵方程时不同？

Question

我一直在处理Python中形式为A = Bx的线性代数问题，并将其与同事在MATLAB和Mathematica中的代码进行比较。当B是一个奇异矩阵时，我们已经注意到Python和其他矩阵之间的区别。当使用numpy.linalg.solve()时，我抛出一个奇异矩阵错误，所以我实现了.pinv()(摩尔彭罗斯伪逆)。

我知道存储逆是计算效率低下的，我首先想知道在Python中是否有更好的方法来处理奇异矩阵。然而，我的问题的重点在于Python如何从无限的解空间中选择一个答案，以及为什么它选择了一个与MATLAB和Mathematica不同的答案。

这是我的玩具问题:

B = np.array([[2,4,6],[1,0,3],[0,7,0]])
A = np.array([[12],[4],[7]])
BI = linalg.pinv(B)
x = BI.dot(A)

Python输出给我的答案是:

[[ 0.4]
[ 1. ]
[ 1.2]]

虽然这肯定是一个正确的答案，但它不是我想要的:(1,1,1)。为什么Python生成这个特殊的解?是否有一种方法可以返回解的空间而不是一个可能的解?我同事的代码返回(1,1,1)- Python不同于Mathematica和MATLAB有什么原因吗?

Answer 1

简而言之，您的代码(显然是np.linalg.lstsq)使用Moore-Penrose伪逆，这在np.linalg.pinv中实现。MATLAB和Mathematica可能会使用高斯消去法来求解系统。我们可以在Python中使用LU分解复制后一种方法:

B = np.array([[2,4,6],[1,0,3],[0,7,0]])
y = np.array([[12],[4],[7]])
P, L, U = scipy.linalg.lu(B)

这将B分解为B = P L U，其中U现在是上对角矩阵，P L是可逆的。特别是，我们发现:

>>> U
array([[ 2.,  4.,  6.],
       [ 0.,  7.,  0.],
       [ 0.,  0.,  0.]])

和

>>> np.linalg.inv(P @ L) @ y
array([[ 12.],
       [  7.],
       [  0.]])

目标是解决这个未确定的转换问题，U x = (P L)^{-1} y。解集和原来的问题是一样的。让解写成x = (x_1, x_2, x_3)。然后我们立即看到，任何解决方案必须有x_2 = 1。那么我们必须有2 x_1 + 4 + 6 x_2 = 12。解x_1，得到x_1 = 4 - 3 x_2。所以任何解的形式都是(4 - 3 x_2, 1, x_2)

生成上述解决方案的最简单方法是简单地选择x_2 = 1。然后是x_1 = 1，然后恢复MATLAB给出的解:(1,1,1)。

另一方面，np.linalg.pinv计算Moore-Penrose伪逆，即满足B伪逆性质的唯一矩阵。这里的重点是唯一的。因此，当你说:

我的问题是Python如何从无限解空间中选择一个答案

答案是，当你使用伪逆时，所有的选择实际上都是由你完成的，因为np.linalg.pinv(B)是一个唯一矩阵，因此np.linalg.pinv(B) @ y是唯一的。

要生成完整的解决方案集，请参阅@ali_m上面的注释