python - numpy.linalg.solve 和 numpy.linalg.lu_solve 之间的区别

Question

要求解线性矩阵方程，可以使用 numpy.linalg.solve它实现了 LAPACK 例程 *gesv .

根据文档

DGESV computes the solution to a real system of linear equations
    A * X = B,
 where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

 The LU decomposition with partial pivoting and row interchanges is
 used to factor A as
    A = P * L * U,
 where P is a permutation matrix, L is unit lower triangular, and U is
 upper triangular.  The factored form of A is then used to solve the
 system of equations A * X = B.

但是，我们也可以使用 scipy.linalg.lu_factor() 和 scipy.linalg.lu_solve()为了解决我们的问题，其中lu_factor()是

Compute pivoted LU decomposition of a matrix.

The decomposition is:

A = P L U

where P is a permutation matrix, 
L lower triangular with unit diagonal elements, 
and U upper triangular.

显然这两种方法似乎是多余的。这种冗余有什么意义吗？似乎让我感到困惑..

Answer 1

确实你是对的:链接 scipy 的 scipy.linalg.lu_factor() 和 scipy.linalg.lu_solve() 完全等同于 numpy 的 numpy.linalg.solve()。不过，在实际情况下访问 LU 分解是一个很大的优势。

首先，让我们证明等价性。 numpy.linalg.solve() 指出:

The solutions are computed using LAPACK routine _gesv

确实，github repository of numpy 包含 LAPACK 的精简版。那么，我们来看看LAPACK的dgesv的源码。计算矩阵的 LU 分解并用于求解线性系统。确实，该函数的源代码非常清晰:用于输入检查的appart，归结为调用dgetrf(LU factorization)和dgetrs。最后，scipy.linalg.lu_factor() 和 scipy.linalg.lu_solve() 分别包装了 dgetrf 和 dgetrs，具有类似 getrf, = get_lapack_funcs(('getrf',), (a1,)) 和 getrs, = get_lapack_funcs(('getrs',), (lu, b1))。

正如@Alexander Reynolds 所注意到的，LU 分解可用于计算矩阵的行列式和秩。实际上，关于行列式，numpy.det 调用 LU 因式分解 _getrf! 请参阅 source of numpy.linalg。然而，numpy 使用 SVD 分解计算矩阵的秩。

但是使用 LU 计算排名并不是将接口(interface)公开给 dgetrf 和 dgetrs 的唯一原因。确实有一些常见的情况，调用一次 dgetrf，将 LU 分解保留在内存中并多次调用 dgetrs 是决定性的优势。以 iterative refinement 为例。必须注意的是，计算 LU 分解比使用分解求解线性系统 (N^2) 花费更多的时间 (N^3)。

让我们看一下 Newton-Raphson method 求解耦合非线性方程组 F(x)=0 其中F: R^N->R^N。执行 Newton-Raphson 迭代需要求解矩阵为雅可比矩阵 J 的线性系统:

其中 x_{i+1} 是未知数。 Jacobian J(x_i) 的计算成本通常很高，更不用说系统必须求解的事实了。因此，通常会考虑准牛顿法，其中建立雅可比行列式的近似值。一个直截了当的想法是不在每次迭代时更新雅可比矩阵，只要残差范数减少就一直迭代。 在这样的过程中，dgetrf会每隔一段时间被调用一次，而dgetrs每次迭代都会调用一次。方案见there。让我们看一个例子:让我们尝试从 x_0=2 开始求解 x^2=1。对于牛顿法的 4 次迭代，我们得到 f(x_4)=9.2e-8 和 f(x_5)<1e-13。但如果雅可比矩阵每十次迭代更新一次，则只需对雅可比矩阵进行两次评估即可得到 f(x_12)=5.7e-11 和 f(x_13)=2.2e-14。

还有更有效的策略，包括每次迭代更新一次 LU 分解而不分解任何矩阵。请参阅 Denis 和 Marwil 的 "Direct Secant Updates of Matrix Factorizations" 以及 Brown 等人的 "EXPERIMENTS WITH QUASI-NEWTON METHODS IN SOLVING STIFF ODE SYSTEMS"。等

因此，在 non-linear finite element analysis 中，并不总是在每次牛顿迭代( option of Code_Aster in paragraph 3.10.2MATRICE 或 diana 或 here 或 there 。