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

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

我一直在尝试Python 3中提供的各种稀疏求解器，并比较它们之间以及与Octave和Matlab的性能。我选择了直接方法和迭代方法，下面将详细解释。

为了生成具有带状结构的适当稀疏矩阵，使用具有N = 250，N = 500和N = 1000的正方形网格的有限元解决了泊松问题。这导致矩阵A = N ^ 2xN ^ 2和向量b = N ^ 2x1的尺寸，即最大NxN为一百万。如果有人有兴趣复制我的结果，我已在下面的链接中上载了矩阵A和向量b（它将在30天后过期）Get systems used here。矩阵存储在三元组I，J，V中，即前两列分别是行和列的索引，第三列是与这些索引对应的值。观察到V中有一些值几乎为零，是故意保留的。尽管如此，在Matlab和Python中执行“间谍”矩阵命令后仍保留了带状结构。

为了进行比较，我使用了以下求解器：

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,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求解器在以下位置进行了测试：

1. Windows。
2. Linux。
3. Mac OS。

条件：

• 在系统解决方案之前和之后立即进行计时。即，不考虑读取矩阵的开销。
• 对每个系统进行十次计时，并计算平均值和标准偏差。

硬件：

• Windows和Linux：Dell intel（R）Core（TM）i7-8850H cpu @ 2.6GHz 2.59GHz，32 Gb RAM DDR4。
• Mac OS：Macbook Pro retina 2014年中期的英特尔（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,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,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％）。

解决方法

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

然后，我尝试使用前述的scipy的直接求解器，intel MKL包装器，Tim Davis的UMFPACK再次使用Matlab，Octave和Python来求解Ax = b。 我的第一个惊喜是Matlab和Octave都可以​​使用A \ b来求解系统，但这并不能确定它是直接求解器，因为它会根据矩阵的特征选择最佳的求解器，请参见{{3} }。但是，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'。 / p>

我把分钟数放在这里。解决时间，

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'）会对解决方案性能产生积极影响。 让我知道是否有人会对这个新系统感兴趣，所以我可以上传它。