如何解决LU 分解速度 vs 传统 Ax=b
了解 LU 分解,以及我使用的教科书声称,一旦计算了初始矩阵 L 和 U(假设没有枢轴矩阵 P 是必需的),解决 Ly = b 然后 Ux = y 比 Ax = 更快/更有效(也就是更少的触发器) b 从头开始。但是,当我使用随机 A ⊆ R^(10 x 10) 矩阵和 b ⊆ R^(10 x 1) 在我的系统上运行它时,结果是慢得多(12.49 秒对 6.17 秒)。
这是我的代码:
import numpy as np
from scipy.linalg import lu
from numpy.random import rand
from numpy.linalg import solve
from timeit import timeit as duration
A,b = rand(500,500),rand(500,1)
P,L,U = lu(A)
t_vanilla = duration("solve(A,b)",setup="from __main__ import solve,A,b")
print(t_vanilla)
t_decomps = duration("solve(U,solve(L,b))",U,b")
print(t_decomps)
关于为什么会这样的任何想法?这是我第一次使用 timeit
库,所以我很可能搞砸了哈哈。我在 Mac OSX、python 3.8 上运行它。
另外,我注意到由于某种原因,scipy.solve
WAY 比 numpy.solve
慢(时间多一倍)......任何关于为什么的直觉?
谢谢!!!
解决方法
在运行 solve(U,solve(L,b))
时,您只需将 np.linalg.solve
两个不同的线性系统后面的 LAPACK gesv 例程交给它盲目求解,而没有考虑 L
和 {{ 的三角结构1}}。这是通过为每个 LU 分解分别执行前向 U
和后向 Ly = b
替换来完成的。这就是为什么您的时间似乎显示的时间大约是单个 Ux = y
调用的两倍。
现在根据你的教科书,更一般地说,已经拥有矩阵 solve
的 LU 分解的数值效率是,如果你有多个右侧 {{1},你可以简单地重用它}} 以避免代价高昂且冗余的分解步骤。为此,您只需要先通过前向代换来求解两个三角系统 A
,然后使用后向代换来解决 b
。让我按如下方式计时:
Ly = b
如果您有多个 Ux = y
需要求解,它确实会更有效率。
请注意,我使用了 scipy.linalg.lu_factor,它存储了 LU 分解的紧凑形式,然后是 scipy.linalg.lu_solve,它使用 import numpy as np
from scipy.linalg import lu_factor,lu_solve,lu
A,b = np.random.rand(500,500),np.random.rand(500,1)
P,L,U = lu(A)
lu_A = lu_factor(A)
# LU + forward backward
%timeit np.linalg.solve(A,b)
>>> 1.4 ms ± 10.6 µs per loop (mean ± std. dev. of 7 runs,1000 loops each)
# 2 x (LU + forward backward)
%timeit np.linalg.solve(U,np.linalg.solve(L,b))
>>> 2.39 ms ± 85.4 µs per loop (mean ± std. dev. of 7 runs,100 loops each)
# forward backward
%timeit lu_solve(lu_A,b)
>>> 69.8 µs ± 6.48 µs per loop (mean ± std. dev. of 7 runs,10000 loops each)
给出的分解来求解系统。
另请注意,当多个 b
不能同时使用时,重用 LU 分解很有用。相反,如果它们随时可用,您可以简单地使用水平堆叠的 lu_factor
进行单个 b
调用,并具有等效的性能,因为求解器本身将为每个剩余的 solve
重用为第一个计算的 LU 分解{1}}。
关于 scipy 的性能问题,如果没有关于 numpy/scipy (b
) 的更多信息,很难说。
希望现在更清楚了。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。