事实上,Numpy 和 Scipy 有效地调用 LAPACK 例程来执行numpy.linalg.inv https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html and scipy.linalg.inv https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.inv.html.
要逆一般矩阵,numpy.linalg.inv https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html solves A.x=np.eye((n,n))
。功能inv() https://github.com/numpy/numpy/blob/v1.17.0/numpy/linalg/linalg.py#L486-L552 calls ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
, which calls https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src#L1693 call_@lapack_func@(¶ms);
where params.B
是单位矩阵并且@lapack_func@
是其中之一sgesv, dgesv, cgesv, zgesv
,它们是一般矩阵的线性求解器。
另一方面,scipy.linalg.inv
calls https://github.com/scipy/scipy/blob/v1.4.1/scipy/linalg/basic.py#L911-L983 getri
, 定义为get_lapack_funcs(('getri'),(a1,))
。它对应于DGETRI() http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f_source.htmllapack 函数,旨在使用 LU 分解计算矩阵的逆,计算公式为DGETRF()
. 因此,如果您正在使用DGETRI()
在 Fortran 中,利用scipy.linalg.inv()
在 python 中可能会实现类似的性能和结果。
大多数 Lapack 函数可以使用以下方式调用scipy.linalg.lapack https://docs.scipy.org/doc/scipy/reference/linalg.lapack.html#module-scipy.linalg.lapack。这是一个使用的示例scipy.linalg.cython_lapack.dgetri()
在 cython 模块中:如何为Python编译C扩展,其中C函数使用LAPACK库? https://stackoverflow.com/questions/52106978/how-to-compile-c-extension-for-python-where-c-function-uses-lapack-library/52132359#52132359下面是一个示例代码,在 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
中提供了链接 DGETRF()+DGETRI() 的 Fortran 代码LAPACK 反演例程奇怪地混合了所有变量 https://stackoverflow.com/questions/26475987/lapack-inversion-routine-strangely-mixes-up-all-variables进行一些更改后,运行:
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 库,则可能会出现不同的性能。要调查此问题,请键入如下命令np.__config__.show()
如图所示将 ATLAS/MKL 链接到已安装的 Numpy https://stackoverflow.com/questions/21671040/link-atlas-mkl-to-an-installed-numpy或报告中的命令如何在 NumPy 和 SciPy 中检查 BLAS/LAPACK 链接? https://stackoverflow.com/questions/9000164/how-to-check-blas-lapack-linkage-in-numpy-and-scipy .
为了进一步利用分布式计算,petsc https://www.mcs.anl.gov/petsc/documentation/faq.html#invertmatrix不鼓励对完整矩阵求逆,因为很少需要这样做。书中还写道MatMatSolve(A,B,X)
, where B
and X
可以使用稠密矩阵来做到这一点。此外,这个函数是在python接口中提供的petsc4py https://www.mcs.anl.gov/petsc/petsc4py-current/docs/apiref/index.html作为方法matSolve(self, Mat B, Mat X)
对于对象petsc4py.PETSc.Mat
。不再维护的元素库 https://github.com/elemental/Elemental被列为实现密集矩阵的直接求解器。虽然 Elemental 库支持 python 接口,但它的分支 Hydrogen 不再支持它。
尽管如此,Elemental 页面还是列出了一些与分布式密集线性代数相关的开源项目。ScaLapack 提供了例程PDGETRI()/PZGETRI() http://www.netlib.org/scalapack/explore-html/d3/df3/pdgetri_8f_source.html使用 LU 分解来反转分布式稠密矩阵。这可能会为更快的反转留下一些空间。