有没有办法使用 python 进一步缩短稀疏求解时间？
Is there a way to further improve sparse solution times using python?
我一直在尝试 Python 3 中可用的不同稀疏求解器，并比较它们之间的性能以及与 Octave 和 Matlab 的性能。我选择了直接和迭代两种方法，我将在下面更详细地解释这一点。
为了生成具有带状结构的适当稀疏矩阵，使用具有 N=250、N=500 和 N=1000 平方网格的有限元求解泊松问题。这导致矩阵A=N^2xN^2和向量b=N^2x1的维度，即最大的NxN是一百万。如果有人有兴趣复制我的结果，我已经在下面上传了矩阵 A 和向量 b link（它将在 30 天后过期）Get systems used here。矩阵存储在三元组 I,J,V 中，即前两列分别是行和列的索引，第三列是这些索引对应的值。请注意，V 中有一些值是故意保留的，这些值几乎为零。尽管如此，在 Matlab 和 Python.
中执行“间谍”矩阵命令后，仍保留了带状结构
为了比较，我使用了以下求解器：
Matlab 和 Octave，直接求解器：规范 x=A\b。
Matlab and Octave, pcg solver: The preconditioned conjugated gradient, pcg solver 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 retina mid 2014 intel (R) quad-core(TM) i7 2.2GHz 16 Gb Ram DDR3.
结果：
观察：
- Matlab A\b 是最快的，尽管使用的是较旧的计算机。
- Linux 和 Windows 版本之间存在显着差异。例如，参见 NxN=1e6 处的直接求解器。尽管 Linux 在 windows (WSL) 下 运行ning。
- 一个人在 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。我知道以前的 posts，例如scipy slow sparse matrix solver
How to use Numba to speed up sparse linear system solvers in Python that are provided in scipy.sparse.linalg?
但我认为 none 和这个一样全面，在使用 python 库时突出了操作系统之间的更多问题。
EDIT_1:
我使用评论中建议的 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))
到原文中提供的脚本post。 “.astype(np.float32)”可以省略，并且该系统的性能会稍微变差（大约 10%）。
我会试着自己回答。为了提供答案，我尝试了一个更苛刻的示例，矩阵大小为 (N,N)，大小约为 50 万乘 50 万，对应的向量为 (N,1)。然而，这比问题中提供的稀疏得多（更密集）。这个以 ascii 格式存储的矩阵约为 1.7 Gb，而示例中的矩阵约为 0.25 Gb（尽管其“大小”更大）。在这里查看它的形状，
然后，我尝试再次使用 Matlab、Octave 和 Python 求解 Ax=b，使用前面提到的来自 scipy 的直接求解器、intel MKL 包装器、来自 Tim Davis 的 UMFPACK。
我的第一个惊喜是Matlab和Octave都可以使用A\b求解系统，这不能肯定它是一个直接求解器，因为它根据矩阵的特性选择最好的求解器，见Matlab's x=A\b。但是，python 的 linalg.spsolve、MKL 包装器和 UMFPACK 在 Windows 和 Linux 中抛出 out-of-memory 错误。在 mac 中，linalg.spsolve 正在以某种方式计算出一个解决方案，尽管它的性能非常差，但它从未出现过内存错误。我想知道内存的处理方式是否因 OS 而异。对我来说，似乎 mac 将内存交换到硬盘驱动器而不是从 RAM 中使用它。与 matlab 相比，Python 中 CG 求解器的性能相当差。然而，为了提高 python 中 CG 求解器的性能，如果首先计算 A=0.5(A+A')（如果显然有一个对称系统），则可以在性能上获得巨大的改进。在 Python 中使用预调节器没有帮助。我尝试使用 sp.linalg.spilu 方法和 sp.linalg.LinearOperator 一起计算预条件子，但性能很差。在matlab中，可以使用不完全Cholesky分解
对于 out-of-memory 问题，解决方案是使用 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') 可以对求解性能产生积极影响。
如果有人对这个新系统感兴趣，请告诉我，我可以上传它。
我一直在尝试 Python 3 中可用的不同稀疏求解器，并比较它们之间的性能以及与 Octave 和 Matlab 的性能。我选择了直接和迭代两种方法，我将在下面更详细地解释这一点。
为了生成具有带状结构的适当稀疏矩阵，使用具有 N=250、N=500 和 N=1000 平方网格的有限元求解泊松问题。这导致矩阵A=N^2xN^2和向量b=N^2x1的维度，即最大的NxN是一百万。如果有人有兴趣复制我的结果，我已经在下面上传了矩阵 A 和向量 b link（它将在 30 天后过期）Get systems used here。矩阵存储在三元组 I,J,V 中，即前两列分别是行和列的索引，第三列是这些索引对应的值。请注意，V 中有一些值是故意保留的，这些值几乎为零。尽管如此，在 Matlab 和 Python.
中执行“间谍”矩阵命令后，仍保留了带状结构为了比较，我使用了以下求解器：
Matlab 和 Octave，直接求解器：规范 x=A\b。
Matlab and Octave, pcg solver: The preconditioned conjugated gradient, pcg solver 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 retina mid 2014 intel (R) quad-core(TM) i7 2.2GHz 16 Gb Ram DDR3.
结果：
观察：
- Matlab A\b 是最快的，尽管使用的是较旧的计算机。
- Linux 和 Windows 版本之间存在显着差异。例如，参见 NxN=1e6 处的直接求解器。尽管 Linux 在 windows (WSL) 下 运行ning。
- 一个人在 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。我知道以前的 posts，例如scipy slow sparse matrix solver
EDIT_1: 我使用评论中建议的 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))
到原文中提供的脚本post。 “.astype(np.float32)”可以省略，并且该系统的性能会稍微变差（大约 10%）。
我会试着自己回答。为了提供答案，我尝试了一个更苛刻的示例，矩阵大小为 (N,N)，大小约为 50 万乘 50 万，对应的向量为 (N,1)。然而，这比问题中提供的稀疏得多（更密集）。这个以 ascii 格式存储的矩阵约为 1.7 Gb，而示例中的矩阵约为 0.25 Gb（尽管其“大小”更大）。在这里查看它的形状，
然后，我尝试再次使用 Matlab、Octave 和 Python 求解 Ax=b，使用前面提到的来自 scipy 的直接求解器、intel MKL 包装器、来自 Tim Davis 的 UMFPACK。
我的第一个惊喜是Matlab和Octave都可以使用A\b求解系统，这不能肯定它是一个直接求解器，因为它根据矩阵的特性选择最好的求解器，见Matlab's x=A\b。但是，python 的 linalg.spsolve、MKL 包装器和 UMFPACK 在 Windows 和 Linux 中抛出 out-of-memory 错误。在 mac 中，linalg.spsolve 正在以某种方式计算出一个解决方案，尽管它的性能非常差，但它从未出现过内存错误。我想知道内存的处理方式是否因 OS 而异。对我来说，似乎 mac 将内存交换到硬盘驱动器而不是从 RAM 中使用它。与 matlab 相比，Python 中 CG 求解器的性能相当差。然而，为了提高 python 中 CG 求解器的性能，如果首先计算 A=0.5(A+A')（如果显然有一个对称系统），则可以在性能上获得巨大的改进。在 Python 中使用预调节器没有帮助。我尝试使用 sp.linalg.spilu 方法和 sp.linalg.LinearOperator 一起计算预条件子，但性能很差。在matlab中，可以使用不完全Cholesky分解
对于 out-of-memory 问题，解决方案是使用 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') 可以对求解性能产生积极影响。 如果有人对这个新系统感兴趣，请告诉我，我可以上传它。