LU 分解速度与传统 Ax=b 的对比
Speed of LU decomposition vs traditional Ax=b
了解 LU 分解，我正在使用的教科书声称一旦计算了初始矩阵 L 和 U (假设不需要旋转矩阵 P），求解 Ly = b 的 faster/more 效率更高（又名更少的触发器），然后Ux = y 比 Ax = b 从头开始。但是，当我 运行 在我的系统上使用 运行dom 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))", setup="from __main__ import solve, L, 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例程两个不同的线性系统，它盲目求解而不考虑[=的三角结构13=] 和 U。这是通过分别对每个 LU 分解执行前向 Ly = b 和后向 Ux = y 替换来完成的。这就是为什么您的时间显示大约是单个 solve 调用时间的两倍。
现在根据你的教科书，更一般地说，来自矩阵 LU 分解的数值效率 A 是，如果你有多个右侧 A，你可以简单地重用它 b 以避免昂贵且多余的因式分解步骤。为此，您只需要通过正向替换解决两个三角系统 Ly = b，然后使用反向替换解决 Ux = y。让我按以下时间计时：
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 来解决，它确实更有效率。
请注意，我使用了 scipy.linalg.lu_factor, which stores a compact form of the LU factorization, followed by scipy.linalg.lu_solve，它使用 lu_factor 给出的分解求解系统。
另请注意，当多个 b 不能同时使用时，重用 LU 分解很有用。相反，如果它们很容易获得，您可以简单地使用水平堆叠 b 进行单个 solve 调用，具有等效的性能，因为求解器本身将为每个剩余 b.
关于 scipy 的性能问题，如果没有关于 numpy/scipy (np.show_config()) 的更多信息，很难说。
希望现在更清楚了。
了解 LU 分解，我正在使用的教科书声称一旦计算了初始矩阵 L 和 U (假设不需要旋转矩阵 P），求解 Ly = b 的 faster/more 效率更高（又名更少的触发器），然后Ux = y 比 Ax = b 从头开始。但是，当我 运行 在我的系统上使用 运行dom 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))", setup="from __main__ import solve, L, 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例程两个不同的线性系统，它盲目求解而不考虑[=的三角结构13=] 和 U。这是通过分别对每个 LU 分解执行前向 Ly = b 和后向 Ux = y 替换来完成的。这就是为什么您的时间显示大约是单个 solve 调用时间的两倍。
现在根据你的教科书，更一般地说，来自矩阵 LU 分解的数值效率 A 是，如果你有多个右侧 A，你可以简单地重用它 b 以避免昂贵且多余的因式分解步骤。为此，您只需要通过正向替换解决两个三角系统 Ly = b，然后使用反向替换解决 Ux = y。让我按以下时间计时：
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 来解决，它确实更有效率。
请注意，我使用了 scipy.linalg.lu_factor, which stores a compact form of the LU factorization, followed by scipy.linalg.lu_solve，它使用 lu_factor 给出的分解求解系统。
另请注意，当多个 b 不能同时使用时，重用 LU 分解很有用。相反，如果它们很容易获得，您可以简单地使用水平堆叠 b 进行单个 solve 调用，具有等效的性能，因为求解器本身将为每个剩余 b.
关于 scipy 的性能问题，如果没有关于 numpy/scipy (np.show_config()) 的更多信息，很难说。
希望现在更清楚了。