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

Question

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

• scipy.linalg.lstsq （时间/剩余）：14秒，626.982
  
• scipy.sparse.linalg.lsmr （时间/残差）：4.5s，626.982（相同精度）
  

现在我已经观察到，当我没有秩不足的情况时，形成正态方程并使用cholesky因子分解来求解是解决我的问题的最快方法。所以我的问题是，如果我对最小范数解不感兴趣，那么当A ^ TA为奇数时，是否有办法获得（A ^ TAx = b）的解。我已经尝试过scipy.linalg.solve，但它为奇异矩阵提供了LinAlgError。我也想知道A是否满足m >> n，条件不良，可能不完全排序的问题，那么就时间，残差精度（或任何其他度量）而言，应该使用哪种方法。任何想法和帮助，我们将不胜感激。谢谢！

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)

但是，对于 size 的矩阵(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。希望一切都有意义。