有没有一种方法可以使用python进一步缩短稀疏的解决方案时间？

Question

我一直在尝试Python3中可用的不同稀疏解算器，并将它们之间的性能以及与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中使用“spy”矩阵命令后，带状结构仍然保持不变

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

Matlab和倍频程，直接解算器：规范的x=A\b

Matlab和倍频程，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解算器在以下环境中进行了测试：

1. 窗户
2. Linux
3. 苹果操作系统

条件：

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

硬件：

• 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内存DDR3

结果:

意见：

• Matlab A\b是最快的，尽管它使用的是较旧的计算机
• Linux和Windows版本之间存在显著差异。例如，请参见NxN=1e6处的直接解算器。尽管Linux是在windows（WSL）下运行的，但情况仍然如此
• 在Scipy解算器中可能存在巨大的分散性。也就是说，如果同一个解决方案运行多次，其中一次可能增加两倍以上
• python中最快的选项可能比在更有限的硬件中运行的Matlab慢近四倍。真的吗

如果您想复制测试，我在这里留下非常简单的脚本。 对于matlab/倍频程：

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更好。欢迎提出建议

另外，我知道以前的职位，例如scipy slow sparse matrix solverIssues using the scipy.sparse.linalg linear system solversHow to use Numba to speed up sparse linear system solvers in Python that are provided in scipy.sparse.linalg? 但我认为没有一个比这一个更全面，这突出了在使用python库时操作系统之间的更多问题

编辑_1： 我使用“英特尔MKL”的QR解算器添加了一个带有结果的新绘图，该解算器使用python包装器，如评论中所建议的。然而，这仍然落后于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）”，性能会变得更差（大约t 10%）用于该系统。

Answer 1

我会尽力回答自己的问题。为了给出答案，我尝试了一个更为苛刻的例子，使用一个大小为（N，N）的矩阵，大约50万乘50万，以及相应的向量（N，1）。然而，这比问题中提供的要稀疏得多（更密集）。这个存储在ascii中的矩阵约为1.7GB，而示例中的矩阵约为0.25GB（尽管其“大小”更大）。看看它的形状

然后，我再次尝试使用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.

所以答案是肯定的，有很多方法可以提高scipy的求解时间。如果工作站内存允许，强烈建议使用UMFPACK（Linux）或“英特尔MKL QR解算器”的包装器。否则，如果要处理对称系统，在使用共轭梯度解算器之前执行A=0.5（A+A'）可以对解决方案性能产生积极影响。 如果有人对这个新系统感兴趣，请告诉我，这样我就可以上传它了