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

Question

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

scipy.linalg.lstsq (time/residual) : 14s, 626.982
scipy.sparse.linalg.lsmr (time/residual) : 4.5s, 626.982 (相同精度)

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

Answer 1

我会说 "correct" 解决这个问题的方法是使用 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.[=可能会节省时间33=]

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

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

Fastest way to solve least square for overdetermined system

python

numpy

linear-algebra

scipy

回复评论