python - 解决超定系统最小二乘法的最快方法

Question

我有一个大小为 m*n 的矩阵 A(m 阶约为 100K，n 约为 500)和一个向量 b。此外，我的矩阵是病态的和秩亏的。现在我想找出 Ax = b 的最小二乘解，为此我比较了一些方法:

• scipy.linalg.lstsq(时间/残差):14s，626.982
  
• scipy.sparse.linalg.lsmr(时间/残差):4.5s，626.982(相同精度)
  

现在我观察到，当我没有形成正规方程的秩亏情况时，使用 cholesky 分解求解它是解决我的问题的最快方法。所以我的问题是，如果我对最小范数解不感兴趣，那么当 A^TA 是单数时，有没有办法得到 (A^TAx=b) 的解(任意)。我试过 scipy.linalg.solve 但它给出了奇异矩阵的 LinAlgError 。此外，我想知道 A 是否满足 m>>n、条件不良、可能不是完整的 col-rank，那么在时间、剩余精度(或任何其他指标)方面应该使用哪种方法。非常感谢任何想法和帮助。谢谢!

Answer 1

我想说解决这个问题的“正确”方法是使用 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' 标志实际上可以给您带来很多好处。你当然可以在这里这样做，因为 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 和一个我自己的 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)

我的 a 不是完全单一的，而是接近单一的。

回复评论

我猜你想做的是在一些线性等式约束下找到一个可行点。这里的问题是您不知道哪些约束是重要的。 A 的 100,000 行中的每一行都给你一个新的约束，其中最多 500 个，但可能要少得多(因为不确定性)，实际上很重要。 SVD 为您提供了一种确定哪些维度重要的方法。我不知道还有其他方法可以做到这一点:您可能会在凸优化或线性规划文献中找到一些东西。如果你先验知道 A 的秩是 r ，那么你可以尝试只找到第一个 r 奇异值和对应的向量，如果 r << n 可能会节省时间。

关于您的其他问题，最小范数解决方案不是“最佳”甚至“正确”解决方案。由于您的系统未定，您需要加入一些额外的约束或假设，这将帮助您找到唯一的解决方案。最小范数约束就是其中之一。最小范数解决方案通常被认为是“好的”，因为如果 x 是您尝试设计的某个物理信号，那么具有较低范数的 x 通常对应于具有较低能量的物理信号，然后转化为成本节约等。或者，如果 x 是您尝试估计的某个系统的参数，那么选择最小范数解决方案意味着您假设系统在某种程度上是高效的，并且只使用产生结果 b 所需的最少能量。希望一切都有意义。