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

Question

我一直在用 Python 处理形式为 A = Bx 的线性代数问题，并将其与同事在 MATLAB 和 Mathematica 中的代码进行比较。当 B 是奇异矩阵时，我们注意到 Python 与其他语言之间的差异。使用 numpy.linalg.solve() 时，我抛出了奇异矩阵错误，因此我改为实现了 .pinv()(Moore Penrose 伪逆)。

我知道存储逆矩阵在计算上效率很低，首先我很好奇在 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 可能使用高斯消元来求解系统。我们可以使用 LU 分解在 Python 中复制后一种方法:

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 的伪逆属性的唯一矩阵.这里的重点是独特。因此，当你说:

my question lies in how Python chooses an answer from an infinite solution space

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

要生成全套解决方案，请参阅@ali_m 上面的评论。