LU分解的必要性(以numpy为例)

Question

我试图理解使用numpy和and库进行LU分解的必要性。据我所知，我们想要解Ax = b，首先将A分解成两个三角矩阵L和U，然后通过求解Ly =b来求解LUx =b，然后求解Ux = y，通过求解三角形矩阵，我们可以减少与高斯消除相比的时间。

所以，我在蟒蛇中厌倦了这个想法，用的是numpy和scipy。

我首先用玩具例子构造A和b：

A = np.array([[2, 1, 0, 5], [1, 2, 1, 2], [0, 1, 2, 4], [1, 3, 6, 4.5]])
b = np.array([9, 10, -2, 3])

然后首先在np.solve中解决这个玩具示例

%timeit np.linalg.solve(A, b )

时间是

每环9.76s±782 ns (平均±std )。dev.7次运行中，每一次循环100000次)

然后，我使用因式分解来解决这个系统：

lu, piv = linalg.lu_factor(A)
%timeit linalg.lu_solve((lu, piv), b)

我看到输出是

18.8±213 ns /环(平均±std )。dev.7次运行中，每一次循环100000次)

，与np.solve相比，这是比较安静的慢。

所以，我的问题是，为什么np.solve比linalg.lu_factor更快？我猜numpy.solve没有用高斯消去来解方程吗？和这里的结果有点混淆。

编辑

现在，我用一个更大的矩阵来做实验(10000 X 10000)。

结果如下:对于np.linalg.solve

8.64 s ± 180 ms per loop (mean ± std. dev. of 7 runs, 1 loop each);

对于scipy.linalg.lu_solve

121 ms ± 3.79 ms per loop (mean ± std. dev. of 7 runs, 10 loops each).

对于lu_solve，我只计算求解时间，不计算分解部分。现在更快了！

Answer 1

这是一个部分的答案，因为我质疑你的一个前提。

你写到"LU解决应该比高斯消除更快“。你似乎误解了LU分解的目的。如果你只解决一个这样的问题(给出矩阵Ax=b和向量b )，LU分解不会比高斯消去更快。事实上，分解的算法非常类似于消除，而且也不会更快。

LU分解的优点在于，当你得到矩阵A时，你想要求解多个不同的给定向量b的方程Ax=b。高斯消除需要从头开始，每个解决方案将花费相同的时间。在LU分解中，您可以从第一次计算中存储得到的矩阵L和U，这大大加快了使用不同向量b的后续方程的求解速度。

您可以在C中关于LU分解及其应用的数值计算部分中更多地阅读这方面的内容。

Answer 2

查看numpy.linalg.solve的docstring。它在"Notes“部分中说，”解决方案是使用LAPACK例程_gesv计算的“。(下划线是对应于数据类型的字符的位置保持符。例如，dgesv使用双重精度。)

dgesv解释说它使用了LU分解。因此，您或多或少地在复制计算，但是在Python中执行更多的步骤，所以代码要慢一些。