在Python中求解表面上（但实际上不是）超定的稀疏线性系统

Question

我有一个稀疏矩阵 A（使用 scipy.sparse），还有一个向量 b，我想解方程 Ax = b，求出 x。矩阵 A 的行数比列数多，看起来像是一个超定方程；不过，A 的行之间是线性相关的，这意味着实际上 A 的行秩等于列数。举个例子，A 可能是这样的：

A = np.array([[1., 1.], [-1., -1.], [1., 0.]])

而 b 是：

b = np.array([0., 0., 1.])

最终的解是 x = [1., -1.]。我想知道如何在 Python 中使用 scipy.sparse.linalg 提供的函数来解决这个系统。谢谢！

Answer 1

解决你问题的正式正确方法是使用奇异值分解（SVD）。你有一个这样的系统：

A [MxN] * x [Nx1] = b [Mx1]

奇异值分解会把矩阵 A 分解成三个其他的矩阵，所以你会得到：

U [MxM] * S[MxN] * V[N*N] * x[Nx1] = b[Mx1]

矩阵 U 和 V 都是正交的（它们的逆就是它们的转置），而 S 是一个对角矩阵。如果我们把上面的内容重新写一下，我们得到：

S[MxN] * V [N * N] * x[Nx1] = U.T [MxM] * b [Mx1]

如果 M > N，那么矩阵 S 的最后 M - N 行会全是零。如果你的系统确实是确定的，那么 U.T b 的最后 M - N 行也应该是零。这意味着你可以这样来解决你的系统：

>>> a = np.array([[1., 1.], [-1., -1.], [1., 0.]])
>>> b = np.array([0., 0., 1.])
>>> u, s, v = np.linalg.svd(a)
>>> np.allclose(u.T.dot(b)[-m+n:], 0) #check system is not overdetermined
True
>>> np.linalg.solve(s[:, None] * v, u.T.dot(b)[:n])
array([ 1., -1.])

Answer 2

你的系统可能是欠定的吗？如果不是，并且确实有解，那么最小二乘解就是那个解，所以你可以尝试一下

from scipy.sparse.linalg import lsqr
return_values = lsqr(A, b)
x = return_values[0]

如果你的系统确实是欠定的，这个方法应该能找到最小L2范数的解。如果不行，可以把参数damp设置得非常小（比如1e-5）。

如果你的系统是完全确定的（也就是说A的秩是满的）并且有解，而且你的矩阵A是高的，正如你所描述的，那么你可以在正常方程中找到一个等效的系统：

A.T.dot(A).dot(x) == A.T.dot(b)

这个系统在x上有唯一解。这是一个方形线性系统，因此可以使用线性系统求解器，比如scipy.sparse.linalg.spsolve来解决。