问题描述
了解 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) 的更多信息，很难说。
希望现在更清楚了。