【问题标题】：The necessity of LU decomposition (using numpy as an example)LU分解的必要性（以numpy为例）
【发布时间】：2018-09-17 10:40:59
【问题描述】：

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

所以，我在 python 中使用 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.76 µs ± 782 ns（7 次运行的平均值 ± 标准偏差，每次 100000 次循环）

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

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

我看到输出是

每个循环 18.8 µs ± 213 ns（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，我只计算求解的时间，分解部分不计算。现在速度更快了！

【问题讨论】：

我想知道“通过求解三角矩阵，与高斯消元法相比，我们可以减少时间。”这一说法的来源。 AFAIU LU分解基本上是高斯消除的一种修改形式，因此LU分解本身的时间复杂度不能比高斯消除更好。如果你想用不同的bs 解决许多Ax = b，LU 分解可能是有益的，但对于单个系统，我认为没有理由让它更快。
尝试更多、很多个更大的矩阵，以确保 Python 与 C 算法对话的开销（相对）尽可能小。在 4x4 矩阵上进行比较毫无意义，为函数调用分配所有帧堆栈将比实际计算花费更长的时间。
@SergGr 也许我之前的问题令人困惑。我是说在LU分解之后，使用L和U来求解方程应该比高斯消去更快。
@AndreyTyukin 我现在使用更大的矩阵来执行此操作，lu_solve 的结果要快得多。

【解决方案1】：

这是一个部分答案，因为我对你的一个前提提出异议。

您写道“LU 求解应该比高斯消元法更快”。您似乎误解了 LU 分解的目的。如果您只解决一个这样的问题（Ax=b，其中给出了矩阵A 和向量b），LU decomp 并不比高斯消除快。确实，分解的算法与消除算法非常相似，而且速度并不快。

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

您可以在 C 中的数字食谱中关于 LU Decomposition and Its Applications 的部分阅读更多相关信息。

【讨论】：

【解决方案2】：

查看numpy.linalg.solve 的文档字符串。它在“注释”部分中说“解决方案是使用 LAPACK 例程 _gesv 计算的”。（下划线是对应于数据类型的字符的占位符。例如，dgesv 使用双精度。）

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

【讨论】：