python - 有没有办法使用python进一步缩短稀疏解决时间？

Question

我一直在尝试使用Python 3中可用的不同稀疏求解器，并比较它们之间以及与Octave和Matlab的性能。我选择了直接方法和迭代方法，下面将详细解释。
为了生成具有带状结构的适当稀疏矩阵，使用具有N = 250，N = 500和N = 1000的正方形网格的有限元解决了泊松问题。这导致矩阵A = N ^ 2xN ^ 2和向量b = N ^ 2x1的尺寸，即最大NxN为一百万。如果有人有兴趣复制我的结果，请在下面的链接Get systems used here中上载矩阵A和向量b(它将在30天后过期)。矩阵存储在三元组I，J，V中，即前两列分别是行和列的索引，第三列是与这些索引对应的值。观察到V中有一些值几乎为零，是故意保留的。尽管如此，在Matlab和Python中使用“ spy ”矩阵命令后仍保留了带状结构。
为了进行比较，我使用了以下求解器:
Matlab和Octave，直接求解器:规范的x=A\b。
Matlab和Octave pcg求解器:预处理的共轭梯度pcg求解器pcg(A,b,1e-5,size(b,1))(不使用预处理器)。
Scipy(Python)，直接求解器:linalg.spsolve(A, b)其中A先前以csr_matrix格式格式化。
Scipy(Python)，PCG解算器:sp.linalg.cg(A, b, x0=None, tol=1e-05)Scipy(Python)，UMFPACK求解器:使用spsolve(A, b)的from scikits.umfpack import spsolve。该求解器显然是在Linux下可用的(仅？)，因为它使用了libsuitesparse [Timothy Davis，Texas A＆M]。在ubuntu中，必须首先将其安装为sudo apt-get install libsuitesparse-dev。
此外，上述python求解器已在以下环境中进行了测试:

Windows。

Linux。

Mac OS。

使适应:

定时在系统解决方案之前和之后完成。即，不考虑读取矩阵的开销。

每个系统的计时完成十次，并计算平均值和标准偏差。

硬件:

Windows和Linux:Dell intel(R)Core(TM)i7-8850H CPU @ 2.6GHz 2.59GHz，32 Gb RAM DDR4。

Mac OS:Macbook Pro视网膜2014年中期intel(R)四核(TM)i7 2.2GHz 16 Gb Ram DDR3。

结果:

观察结果:

Matlab A\b是最快的，尽管它使用的是较旧的计算机。

Linux和Windows版本之间存在显着差异。例如，参见NxN = 1e6处的直接求解器。尽管Linux在Windows(WSL)下运行，这仍然存在。

在Scipy求解器中可能会有很大的分散。这就是说，如果同一解决方案运行多次，其中一次可能会增加两倍以上。

python中最快的选项可能比在有限硬件中运行的Matlab慢近四倍。真的吗？

如果您想重现测试，我在这里留下了非常简单的脚本。
对于matlab/octave:

IJS=load('KbN1M.txt');
b=load('FbN1M.txt');

I=IJS(:,1);
J=IJS(:,2);
S=IJS(:,3);

Neval=10;
tsparse=zeros(Neval,1);
tsolve_direct=zeros(Neval,1);
tsolve_sparse=zeros(Neval,1);
tsolve_pcg=zeros(Neval,1);
for i=1:Neval
    tic
    A=sparse(I,J,S);
    tsparse(i)=toc;
    tic
    x=A\b;
    tsolve_direct(i)=toc;        
    tic
    x2=pcg(A,b,1e-5,size(b,1));
    tsolve_pcg(i)=toc;
end

save -ascii octave_n1M_tsparse.txt tsparse
save -ascii octave_n1M_tsolvedirect.txt tsolve_direct
save -ascii octave_n1M_tsolvepcg.txt tsolve_pcg

对于python:

import time
from scipy import sparse as sp
from scipy.sparse import linalg
import numpy as np
from scikits.umfpack import spsolve, splu #NEEDS LINUX


b=np.loadtxt('FbN1M.txt')
triplets=np.loadtxt('KbN1M.txt')

I=triplets[:,0]-1
J=triplets[:,1]-1
V=triplets[:,2]

I=I.astype(int)
J=J.astype(int)
NN=int(b.shape[0])

Neval=10
time_sparse=np.zeros((Neval,1))
time_direct=np.zeros((Neval,1))
time_conj=np.zeros((Neval,1))
time_umfpack=np.zeros((Neval,1))
for i in range(Neval):
    t = time.time()
    A=sp.coo_matrix((V, (I, J)), shape=(NN, NN))
    A=sp.csr_matrix(A)
    time_sparse[i,0]=time.time()-t
    t = time.time()
    x=linalg.spsolve(A, b)
    time_direct[i,0] = time.time() - t
    t = time.time()
    x2=sp.linalg.cg(A, b, x0=None, tol=1e-05)
    time_conj[i,0] = time.time() - t
    t = time.time()
    x3 = spsolve(A, b) #ONLY IN LINUX
    time_umfpack[i,0] = time.time() - t

np.savetxt('pythonlinux_n1M_tsparse.txt',time_sparse,fmt='%.18f')
np.savetxt('pythonlinux_n1M_tsolvedirect.txt',time_direct,fmt='%.18f')
np.savetxt('pythonlinux_n1M_tsolvepcg.txt',time_conj,fmt='%.18f')
np.savetxt('pythonlinux_n1M_tsolveumfpack.txt',time_umfpack,fmt='%.18f')

有没有办法使用python进一步缩短稀疏解决时间？或至少在性能上与Matlab相似？我愿意接受使用C/C++或Fortran和python包装器的建议，但我相信它不会比UMFPACK选择好得多。建议是非常欢迎的。
P.S.我知道以前的帖子，例如 scipy slow sparse matrix solver
Issues using the scipy.sparse.linalg linear system solvers
How to use Numba to speed up sparse linear system solvers in Python that are provided in scipy.sparse.linalg?
但是我认为没有一个比这更全面了，在使用python库时，突出了操作系统之间的更多问题。
EDIT_1:
我使用注释中建议的使用python包装器的intel MKL的QR解算器添加了带有结果的新绘图。但是，这仍然落后于Matlab的性能。
为此，需要添加:

from sparse_dot_mkl import sparse_qr_solve_mkl

和

sparse_qr_solve_mkl(A.astype(np.float32), b.astype(np.float32))

原始帖子中提供的脚本。可以省略“.astype(np.float32)”，并且此系统的性能变差(大约10％)。

Answer 1

我会尽力回答自己的问题。为了提供答案，我尝试了一个要求更高的示例，该示例的大小为(N，N)的矩阵约为50万乘以100万，并具有相应的向量(N，1)。但是，这比问题中提供的要稀疏得多(更密集)。与示例中的矩阵相比，存储在ascii中的此矩阵约为1.7 Gb(尽管其“大小”较大)约为0.25 Gb。在这里看到它的形状，

然后，我尝试再次使用Matlab，Octave和Python，使用前述的scipy的直接求解器，intel MKL包装器，Tim Davis的UMFPACK来求解Ax = b。
我的第一个惊喜是Matlab和Octave都可以​​使用A\b来求解系统，但这并不能确定它是直接求解器，因为它会根据矩阵的特征选择最佳的求解器，请参见Matlab's x=A\b。但是，Python的linalg.spsolve，MKL包装器和UMFPACK在Windows和Linux中抛出了内存不足错误。在Mac中，linalg.spsolve某种程度上在计算解决方案，而性能却很差，它永远不会出现内存错误。我想知道是否根据操作系统对内存进行了不同的处理。在我看来，mac似乎将内存交换到了硬盘驱动器上，而不是从RAM上使用了它。与matlab相比，Python中CG求解器的性能相当差。但是，为了提高python中CG求解器的性能，如果首先计算A = 0.5(A + A')(如果显然有一个对称系统)，则可以大大提高性能。在Python中使用前置条件没有帮助。我尝试将sp.linalg.spilu方法与sp.linalg.LinearOperator一起使用来计算前置条件，但性能相当差。在matlab中，可以使用不完整的Cholesky分解。
对于内存不足问题，解决方案是使用LU分解并求解两个嵌套系统，例如Ax = b，A = LL'，y = L\b和x = y\L'。
我把这里的分钟。解决时间，

Matlab mac, A\b = 294 s.
Matlab mac, PCG (without conditioner)= 17.9 s.
Matlab mac, PCG (with incomplete Cholesky conditioner) = 9.8 s.
Scipy mac, direct = 4797 s.
Octave, A\b = 302 s.
Octave, PCG (without conditioner)= 28.6 s.
Octave, PCG (with incomplete Cholesky conditioner) = 11.4 s.
Scipy, PCG (without A=0.5(A+A'))= 119 s.
Scipy, PCG (with A=0.5(A+A'))= 12.7 s.
Scipy, LU decomposition using UMFPACK (Linux) = 3.7 s total.

因此答案是肯定的，有多种方法可以缩短科学求解时间。如果工作站的内存允许，强烈建议将包装器用于UMFPACK(Linux)或intel MKL QR解算器。否则，如果要处理对称系统，则在使用共轭梯度求解器之前执行A = 0.5(A + A')会对解决方案性能产生积极影响。
让我知道是否有人会对这个新系统感兴趣，所以我可以上传它。