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()指出：

使用LAPACK例程_gesv计算解决方案

实际上，numpy的github存储库包含LAPACK的精简版本。然后，让我们看一下dgesvLAPACK 的来源。计算矩阵的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！请参阅numpy.linalg的来源。但是，numpy使用SVD分解计算矩阵的秩。

但是，使用LU计算等级并不是将接口公开给dgetrfand 的唯一原因dgetrs。实际上，在常见情况下，dgetrf一次校准，将LU因式分解保留在内存中并dgetrs多次调用是决定性的优势。例如，看一下迭代优化。必须注意，计算LU分解比使用线性分解求解线性系统（N ^ 2）要花费更多的时间（N ^ 3）。

让我们看一下牛顿-拉夫森方法来求解非线性方程 F（x）= 0的系统，其中F：R ^ N-> R ^ N。执行牛顿-拉夫森迭代需要求解线性系统，其中矩阵是雅可比矩阵J：

其中x_ {i + 1}是未知的。雅可比定律J（x_i）通常计算起来很昂贵，而不是提及必须解决系统的事实。结果，通常考虑拟牛顿法，其中建立了雅可比方程的近似值。一个简单的想法是，只要残差范数减少，就不必每次迭代都更新Jacobian矩阵并保持迭代。在这样的过程中，dgetrf有时会调用一次，而dgetrs每次迭代都会调用一次。看那里一个计划。让我们看一个例子：让我们尝试从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的“ 矩阵因式分解的直接割线更新”和Brown等人的“用拟牛顿法求解斯蒂夫ODE系统的实验”。等

其结果是，在非线性有限元分析，该组件的切线刚度并不总是评价每牛顿迭代（在段落3.10.2矩阵的计算Code_Aster的选项或黛安娜或这里或那里。