使用 python 求解 7000x7000 线性系统时的最佳性能方法

Question

I\xe2\x80\x99m 需要一种有效的方法来反转 python 中的 7000x7000 空气动力学影响系数（密集）矩阵。在使用 FORTRAN 例程之前，我已经开始使用 LAPACK 中的 LU 分解例程来处理问题，我已经看到它在其他相关应用程序中的使用非常有效。不过，我读到 NumPy 和 SciPy 线性系统求解器大多基于对 C 中相同 LAPACK/BLAS 函数的直接调用，并且想知道切换到 FORTRAN 是否真的会在一定程度上减少计算时间这证明放弃一种更简单、更高级的语言是合理的。

\n\n

如果有 Python 求解器可以保证该大小（1000 到 10000，平方）的矩阵具有相似的性能，那么它们是哪些？

\n\n

我确实需要矩阵逆，所以切换到迭代 Ax=b 解决方案不是一个选择。

\n

Answer 1

事实上，Numpy 和 Scipy 有效地调用 LAPACK 例程来执行numpy.linalg.inv和scipy.linalg.inv。

要逆一般矩阵，numpy.linalg.inv求解A.x=np.eye((n,n))。函数inv()调用ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)，它调用 call_@lapack_func@(&params);其中params.B是单位矩阵 和@lapack_func@是其中之一sgesv, dgesv, cgesv, zgesv，它们是一般矩阵的线性求解器。

另一方面，scipy.linalg.inv 调用 getri，定义为get_lapack_funcs(('getri'),(a1,))。它对应于DGETRI()lapack 函数，旨在使用 计算出的 LU 分解来计算矩阵的逆DGETRF()。因此，如果您在 Fortran 中使用DGETRI()，在 python 中使用scipy.linalg.inv()可能会实现类似的性能和结果。

大多数 Lapack 函数都可以使用 来调用scipy.linalg.lapack。scipy.linalg.cython_lapack.dgetri()以下是在 cython 模块中 使用的示例： How to compile C extension for Python where C function using LAPACK library? 下面是一个示例代码，在 1000x1000 矩阵上比较 scipy.linalg.cython_lapack.dgetrf()+scipy.linalg.cython_lapack.dgetri() 、 numpy 和 scipy.linalg.inv() ：

import numpy as np
from scipy import linalg
import time

import myinverse
n=1000
A=np.random.rand(n,n)

start= time.time()
Am,info,string=myinverse.invert(A.copy())
end= time.time()
print 'DGETRF+DGETRI, ', end-start, ' seconds'
if info==0:
    print 'residual ',np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)
else :
    print "inversion failed, info=",info, string

start= time.time()
Am=np.linalg.inv(A.copy())
end= time.time()
print 'np.linalg.inv ', end-start, ' seconds'
print 'residual ', np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)

start= time.time()
Am=linalg.inv(A.copy())
end= time.time()
print 'scipy.linalg.inv ', end-start, ' seconds'
print 'residual ',np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)

输出是：

DGETRF+DGETRI,  0.22541308403  seconds
residual  4.2155882951089296e-11
np.linalg.inv  0.29932808876  seconds
residual  4.371813154546711e-11
scipy.linalg.inv  0.298856973648  seconds
residual  9.110997546690758e-11

对于 2000x2000 矩阵：

DGETRF+DGETRI,  1.64830899239  seconds
residual  8.541625644634121e-10
np.linalg.inv  2.02795410156  seconds
residual  7.448244269611659e-10
scipy.linalg.inv  1.61937093735  seconds
residual  1.6453560233026243e-09

LAPACK 反演例程中提供了链接 DGETRF()+DGETRI() 的 Fortran 代码，奇怪地混合了所有变量 经过一些更改后，让我们运行：

PROGRAM solvelinear
implicit none
REAL(8), dimension(1000,1000)     :: A,Ainv,M,LU
REAL(8),allocatable              :: work(:)
REAL(8)                          :: wwork
INTEGER                        :: info,lwork
INTEGER,dimension(1000)        :: ipiv
INTEGER                        :: i,j
real :: start, finish

        ! put code to test here


info=0
!work=0
ipiv=0

call RANDOM_NUMBER(A)

call cpu_time(start)
!-- LU factorisation
LU = A
CALL DGETRF(1000,1000,LU,1000,ipiv,info)

!-- Inversion of matrix A using the LU
Ainv=LU
lwork=-1
CALL DGETRI(1000,Ainv,1000,Ipiv,wwork,lwork,info)
lwork =INT( wwork+0.1)
allocate(work(lwork))
CALL DGETRI(1000,Ainv,1000,Ipiv,work,lwork,info)
deallocate(work)

call cpu_time(finish)
print '("Time = ",f6.3," seconds.")',finish-start

!-- computation of A^-1 * A to check the inverse
M = matmul(Ainv,A)

print*,"M = "
do i=1,3
  do j=1,3
    print*,M(i,j)
  enddo
end do

END PROGRAM solvelinear

使用 编译后gfortran main2.f90 -o main2 -llapack -lblas -lm -Wall，1000x1000 矩阵需要 0.42s，2000x2000 矩阵需要 3s。

最后，如果 Fortran 代码和 python 代码不链接到相同的 Blas/Lapack 库，则可能会出现不同的性能。要调查此问题，请键入将ATLAS/MKL 链接到已安装的 Numpynp.__config__.show()中所示的命令或如何检查 NumPy 和 SciPy 中的 BLAS/LAPACK 链接？中报告的命令。。

为了进一步利用分布式计算，petsc不鼓励反转完整矩阵，因为很少需要这样做。还指出MatMatSolve(A,B,X)、 其中B和X是稠密矩阵可以用来做到这一点。此外，这个函数在 python 接口petsc4pymatSolve(self, Mat B, Mat X)中作为对象的方法提供petsc4py.PETSc.Mat。不再维护的Elemental 库被列为实现密集矩阵的直接求解器。虽然 Elemental 库支持 python 接口，但它的分支 Hydrogen 不再支持它。尽管如此，Elemental 页面还是列出了一些与分布式密集线性代数相关的开源项目。ScaLapack 提供了例程PDGETRI()/PZGETRI()使用 LU 分解来反转分布式稠密矩阵。这可能会为更快的反转留下一些空间。