Difference between numpy.linalg.solve and numpy.linalg.lu_solve

Question

To solve a linear matrix equation one can use numpy.linalg.solve which implements LAPACK routine *gesv.

According to the documentation

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.

However, we can also use scipy.linalg.lu_factor() and scipy.linalg.lu_solve() in order to solve our problem, where lu_factor() is

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.

So apparently the two methods seem to be redundant. Is there any sense for this redundancy? Seems confusing to me..

Answer 1

The time complexity of solve is O(n³), while that of lu_solve is O(n²). But that of lu_factor is O(n³). So, when you have many systems Ax=b with the same A but different bs, it is more efficient to lu_factor once and lu_solve than to use solve.

Answer 2

Indeed you are right: chaining scipy's scipy.linalg.lu_factor() and scipy.linalg.lu_solve() is perfectly equivalent to numpy's numpy.linalg.solve(). Nevertheless, having access to the LU decomposition is a great advantage in practical situations.

First, let's proove the equivalence. numpy.linalg.solve() states that:

The solutions are computed using LAPACK routine _gesv

Indeed, the github repository of numpy contains a lite version of LAPACK. Then, let's take a look at the source of dgesv of LAPACK. The LU factorization of the matrix is computed and used to solve the linear system. Indeed, the source code of the function is very clear: appart for input checking, it boils down to calling dgetrf (LU factorization) and dgetrs. Finally, scipy.linalg.lu_factor() and scipy.linalg.lu_solve() respectively wraps dgetrf and dgetrs, featuring lines like getrf, = get_lapack_funcs(('getrf',), (a1,)) and getrs, = get_lapack_funcs(('getrs',), (lu, b1)).

As noticed by @Alexander Reynolds , the LU decomposition can be useful to compute the determinant and the rank of the matrix. Indeed, regarding the determinant, numpy.det call the LU factorization _getrf ! See the source of numpy.linalg. Nevertheless, numpy compute the rank of matrices using the SVD decomposition.

But computing the rank using LU is not the only reason for exposing interfaces to dgetrf and dgetrs. There are indeed common situations where caling dgetrf once, keeping the LU factorization in memory and calling dgetrs many times is a decisive advantage. Look at iterative refinement for instance. It must be noticed that computing the LU factorization takes much more time (N^3) than solving a linear system using the factorization (N^2).

Let's take a look at the Newton-Raphson method to solve a system of coupled non-linear equation F(x)=0 where F: R^N->R^N. Performing a Newton-Raphson iteration requires solving a linear system where the matrix is the Jacobian matrix J:

where x_{i+1} is the unknown. The Jacobian J(x_i) is often expensive to compute, not to mension the fact that the system must be solved. As a result, quasi Newton methods are often considered, where approximations of the Jacobian are built. A straightforward idea is not to update the Jacobian matrix every iteration and keep iterating as long as a norm of the residu decreases. In such a process, dgetrf is called from time to time, while dgetrs is called once for every iteration. See there for a scheme. Let's look at an exemple : let's try to solve x^2=1 starting at x_0=2. For 4 iterations of the Newton method, we get f(x_4)=9.2e-8 and f(x_5)<1e-13. But if the Jacobian is updated once every ten iterations, only two evaluations of the Jacobian are required to get f(x_12)=5.7e-11 and f(x_13)=2.2e-14.

There are even more efficient strategies which consist in updating the LU decomposition once every iteration without factoring any matrix. See "Direct Secant Updates of Matrix Factorizations" of Denis and Marwil and "EXPERIMENTS WITH QUASI-NEWTON METHODS IN SOLVING STIFF ODE SYSTEMS" by Brown et. al.

As a result, in non-linear finite element analysis, the assembly of the tangent stiffness is not always evaluated every Newton iteration (option of Code_Aster in paragraph 3.10.2 MATRICE or diana or here or there.