python-3.x - LU 分解速度与传统 Ax=b 的比较

Question

学习 LU 分解，我使用的教科书声称一旦计算了初始矩阵 L 和 U(假设没有主元矩阵 P 是需要的)，解决Ly = b然后Ux = y比Ax =更快/更高效(也称为更少的失败) b从头开始。但是，当我使用随机 A ⊆ R^(10 x 10) 矩阵和 b ⊆ R^(10 x 1) 在系统上运行它时，结果是慢得多(12.49 秒 vs 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 比 numpy.solve 慢方式(时间超过两倍)。 .有什么直觉可以解释为什么吗？

谢谢!!!

Answer 1

当运行solve(U,solve(L,b))时，你只需给LAPACK gesv np.linalg.solve 背后的例程是盲目求解的两个不同的线性系统，而不考虑 L 和 U 的三角形结构。这是通过分别对每个 LU 分解然后进行前向 (Ly = b) 和后向 (Ux = y) 替换来完成的。这就是为什么您的计时似乎表明大约是单个 solve 调用时间的两倍。

现在，根据您的教科书，更一般地说，已经对矩阵 A 进行 LU 分解所带来的数值效率是，如果您有多个右侧，您可以简单地重用它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 ，它存储 LU 分解的紧凑形式，后跟 scipy.linalg.lu_solve ，它通过 lu_factor 给出的分解来求解系统。

另请注意，当这些多个 b 不能同时可用时，重用 LU 分解非常有用。相反，如果它们很容易获得，您可以简单地使用水平堆叠的 b 进行单个 solve 调用，具有相同的性能，因为求解器本身将重用为剩下的每个 b 的第一个。

关于 scipy 的性能问题，如果没有有关 numpy/scipy 的更多信息(np.show_config())，很难说。

希望现在更清楚了。