解决超定系统的最小平方的最快方法 [英] Fastest way to solve least square for overdetermined system

问题描述

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

I have a matrix A of size m*n( m order of ~100K and n ~500) and a vector b. Also, my matrix is ill-conditioned and rank-deficient. Now I want to find out the least-square solution to Ax = b and to this end I have compared some of the methods:

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

• scipy.linalg.lstsq (time/residual) : 14s, 626.982
  
• scipy.sparse.linalg.lsmr (time/residual) : 4.5s, 626.982 (same accuracy)
  

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

Now I have observed that when I don't have the rank-deficient case forming the normal equation and solving it using cholesky factorization is fastest way to solve my problem. So my question is this If I am not interested in the minimum norm solution then is there a way to get a solution(any) to (A^TAx=b) when A^TA is singular. I have tried scipy.linalg.solve but it gives LinAlgError for singular matrices. Also I would like to know if A is such that m>>n, ill-conditonied, possibly not full col-rank then which method should one use in terms of time, residual accuracy(or any other metric). Any thoughts and help is greatly appreciated. Thanks!

推荐答案

我要说的正确"方法是使用SVD，查看您的奇异值谱，并找出多少奇异值您要保持，即弄清楚您希望A^T x与b的距离.遵循以下原则:

I'd say the "correct" way of going about this is to use the SVD, look at your singular value spectrum, and figure out how many singular values you want to keep, i.e., figure out how close you want A^T x to be to b. Something along these lines:

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)的矩阵，这太慢了.我建议您自己实现最小二乘，并添加少量正则化，以避免矩阵变得奇异的问题.

However, for a matrix of size (100000, 500), this is just going to be way too slow. I would recommend implementing least squares by yourself, and adding a small amount of regularization to avoid the issue of the matrix becoming singular.

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')

这是我的工作站上的时间分析*:

Here's a timing analysis on my workstation*:

>>> %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，所以我无法与之比较.

*I somehow don't seem to have scipy.sparse.linalg.lsmr on my machine, so I couldn't compare against that.

在这里似乎并没有做很多事情，但是我在其他地方看到过添加assume_a='pos'标志实际上可以为您带来很多好处.您肯定可以在这里执行此操作，因为保证A^T A为正半定数，并且lamda使其为正定数.您可能需要稍微玩lamda才能将错误降低到足够低的水平.

It doesn't seem to do much here, but I've seen elsewhere that adding the assume_a='pos' flag can actually give you a lot of benefit. You can certainly do this here, since A^T A is guaranteed to be positive semi-definite, and the lamda makes it positive definite. You might have to play with lamda a little to bring your error sufficiently low.

关于错误:

>>> 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:

PS: I generated an a and a b of my own like this:

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不是完全单数，而是几乎单数.

My a isn't completely singular, but near-singular.

我猜您想做的是在某些线性相等约束下找到一个可行的点.这里的问题是您不知道哪些约束很重要. A的100,000行中的每一行都给您一个新的约束，实际上是最重要的约束，但其中最多500约束，但实际上可能要少得多(由于不确定性). SVD为您提供了确定哪些尺寸重要的方法.我不知道这样做的另一种方式:您可能会在凸优化或线性编程文献中找到一些东西.如果您先验地知道A的等级是r，那么您可以尝试仅查找第一个r奇异值和相应的向量，如果r << n可能会节省时间.

I guess what you are trying to do is to find a feasible point under some linear equality constraints. The trouble here is that you don't know which constraints are important. Each of the 100,000 rows of A gives you a new constraint, out of which at most 500, but possibly far fewer (because of underdetermined-ness), actually matter. The SVD gives you a way of figuring out which dimensions are important. I don't know of another way to do this: you might find something in convex optimization or linear programming literature. If you know a priori that the rank of A is r, then you can try to find only the first r singular values and corresponding vectors, which might save time if r << n.

关于您的其他问题，最低标准解决方案不是最佳"解决方案，甚至不是正确"解决方案.由于系统尚未确定，因此您需要添加一些其他约束或假设，这将有助于您找到独特的解决方案.最小规范约束就是这样一种.最小范数解通常被认为是好的"，因为如果x是您要设计的某种物理信号，则范数较低的x通常对应于能量较低的物理信号，然后转换为节省成本等.或者，如果x是您要估算的某些系统的参数，则选择最小范数解意味着您假设系统在某种程度上是有效的，并且仅使用产生结果b所需的最小能量.希望一切都有道理.

Regarding your other question, the minimum norm solution isn't the "best" or even the "correct" solution. Since your system is underdetermined, you need to throw in some additional constraints or assumptions which will help you find a unique solution. The minimum norm constraint is one such. The minimum norm solution is often considered to be "good", because if x is some physical signal which you are trying to design, then an x with lower norm often corresponds to a physical signal with lower energy, which then translates to cost savings, etc. Alternatively, if x is a parameter of some system you are trying to estimate, then choosing the minimum norm solution means you are assuming that the system is efficient in some way, and uses only the minimum energy needed to produce the outcome b. Hope that all makes sense.

这篇关于解决超定系统的最小平方的最快方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！