我想说解决这个问题的“正确”方法是使用 SVD，查看您的奇异值谱，并找出您想要保留多少个奇异值，即找出您想要的接近程度A^T x成为b。沿着这些思路：
def svd_solve(a, b):
[U, s, Vt] = la.svd(a, full_matrices=False)
r = max(np.where(s >= 1e-12)[0])
temp = np.dot(U[:, :r].T, b) / s[:r]
return np.dot(Vt[:r, :].T, temp)
然而，对于大小的矩阵(100000, 500)，这实在是太慢了。我建议您自己实现最小二乘法，并添加少量正则化以避免矩阵变得奇异的问题。
def naive_solve(a, b, lamda):
return la.solve(np.dot(a.T, a) + lamda * np.identity(a.shape[1]),
np.dot(a.T, b))
def pos_solve(a, b, lamda):
return la.solve(np.dot(a.T, a) + lamda * np.identity(a.shape[1]),
np.dot(a.T, b), assume_a='pos')
这是我的工作站上的时序分析*：
>>> %timeit la.lstsq(a, b)
1.84 s ± 39.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>>> %timeit naive_solve(a, b, 1e-25)
140 ms ± 4.15 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
>>> %timeit pos_solve(a, b, 1e-25)
135 ms ± 768 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
*我似乎没有scipy.sparse.linalg.lsmr在我的机器上，所以我无法与之比较。
它在这里似乎没有多大作用，但我在其他地方看到添加assume_a='pos'flag实际上可以给你带来很多好处。您当然可以在这里执行此操作，因为A^T A保证是半正定的，并且lamda使其正定。你可能需要玩lamda一点点使你的错误足够低。
就错误而言：
>>> xhat_lstsq = la.lstsq(a, b)[0]
>>> la.norm(np.dot(a, xhat_lstsq) - b)
1.4628232073579952e-13
>>> xhat_naive = naive_solve(a, b, 1e-25)
>>> la.norm(np.dot(a, xhat_naive) - b)
7.474566255470176e-13
>>> xhat_pos = pos_solve(a, b, 1e-25)
>>> la.norm(np.dot(a, xhat_pos) - b)
7.476075564322223e-13
PS：我生成了一个a and a b我自己的像这样：
s = np.logspace(1, -20, 500)
u = np.random.randn(100000, 500)
u /= la.norm(u, axis=0)[np.newaxis, :]
a = np.dot(u, np.diag(s))
x = np.random.randn(500)
b = np.dot(a, x)
My a不是完全奇异的，而是接近奇异的。
对评论的回应
我想你想做的是在一些线性等式约束下找到一个可行点。这里的问题是你不知道哪些约束是重要的。 A 的 100,000 行中的每一行都会给您一个新的约束，其中最多 500 个（但可能少得多（因为不确定性））实际上很重要。 SVD 为您提供了一种确定哪些维度重要的方法。我不知道还有其他方法可以做到这一点：您可能会在凸优化或线性编程文献中找到一些东西。如果您先验知道 A 的等级是r，那么你可以尝试只找到第一个r奇异值和相应的向量，如果r << n.
关于您的其他问题，最低范数解决方案不是“最佳”甚至“正确”的解决方案。由于您的系统是不确定的，因此您需要引入一些额外的约束或假设，这将帮助您找到独特的解决方案。最小范数约束就是其中之一。最小范数解通常被认为是“好的”，因为如果x是您尝试设计的一些物理信号，然后是x较低范数通常对应于具有较低的能量，然后转化为成本节省等。或者，如果x是您尝试估计的某个系统的参数，那么选择最小范数解意味着您assuming该系统在某种程度上是高效的，并且仅使用产生结果所需的最少能量b。希望一切都有意义。
