有没有办法使用 python 进一步改善稀疏求解时间？

Question

我一直在尝试 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, 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. Mac 操作系统。

条件：

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

硬件：

• Windows 和 Linux：Dell intel (R) Core(TM) i7-8850H CPU @2.6GHz 2.59GHz，32 Gb RAM DDR4。
• Mac 操作系统：Macbook Pro Retina 2014 年中 intel (R) quad-core(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 选择更好。非常欢迎提出建议。

附：我知道以前的帖子，例如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： 我使用 cmets 中建议的 python 包装器使用来自英特尔 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 万乘以 50 万，并且对应的向量为 (N,1)。但是，这比问题中提供的要少得多（更密集）。存储在 ascii 中的这个矩阵大约为 1.7 Gb，而示例中的矩阵大约为 0.25 Gb（尽管它的“大小”更大）。在这里看看它的形状，

然后，我尝试使用上述来自 scipy 的直接求解器、英特尔 MKL 包装器、来自 Tim Davis 的 UMFPACK 再次使用 Matlab、Octave 和 Python 求解 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) 或 intel MKL QR 求解器的包装器。否则，如果处理对称系统，在使用共轭梯度求解器之前执行 A=0.5(A+A') 会对求解性能产生积极影响。 如果有人有兴趣拥有这个新系统，请告诉我，以便我上传。