python – 解决超定系统最小二乘法的最快方法
内容导读
互联网集市收集整理的这篇技术教程文章主要介绍了python – 解决超定系统最小二乘法的最快方法，小编现在分享给大家，供广大互联网技能从业者学习和参考。文章包含3308字，纯文字阅读大概需要5分钟。
内容图文
我有一个大小为m * n(m级为~100K且n~500)的矩阵A和一个向量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,那么在时间,剩余准确度(或任何其他度量)方面应该使用哪种方法.非常感谢任何想法和帮助.谢谢！
解决方法:
我要说的是“正确”的方法是使用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 <<< ñ. 关于您的其他问题,最小规范解决方案不是“最佳”甚至是“正确”的解决方案.由于您的系统不确定,您需要引入一些额外的约束或假设,这将有助于您找到一个独特的解决方案.最小范数约束就是这样.最小范数解决方案通常被认为是“好的”,因为如果x是您尝试设计的某些物理信号,那么具有较低范数的x通常对应于具有较低能量的物理信号,然后转化为成本节省,或者,如果x是您尝试估算的某个系统的参数,那么选择最小范数解决方案意味着您假设系统在某种程度上是有效的,并且仅使用产生结果b所需的最小能量.希望一切都有道理.
内容总结
以上是互联网集市为您收集整理的python – 解决超定系统最小二乘法的最快方法全部内容，希望文章能够帮你解决python – 解决超定系统最小二乘法的最快方法所遇到的程序开发问题。 如果觉得互联网集市技术教程内容还不错，欢迎将互联网集市网站推荐给程序员好友。
内容备注
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 gblab@vip.qq.com 举报，一经查实，本站将立刻删除。
内容手机端
扫描二维码推送至手机访问。