计算科学 - 为什么稀疏矩阵的“scipy.sparse.linalg.spilu”效率低于“scipy.linalg.lu”？ - 吾爱随笔录

为什么稀疏矩阵的“scipy.sparse.linalg.spilu”效率低于“scipy.linalg.lu”？

计算科学 Python 稀疏矩阵 scipy 矩阵分解

2021-12-16 01:11:24

我有一个B稀疏矩阵，并尝试使用scipy.sparse.linalg.spilu专门用于稀疏矩阵的函数来分解B。你能解释一下为什么这个函数比scipy.linalg.lu一般矩阵的函数效率低吗？太感谢了！

import numpy as np
import scipy.linalg as la
import scipy.sparse.linalg as spla
import time
from scipy import sparse
from scipy.sparse import csc_matrix
A = np.random.randint(100, size=(10000, 10000))
B = np.triu(A, -100)

start = time.time()
(P, L, U) = la.lu(B)
end = time.time()
print('Time to decompose B with lu =', end - start)

start = time.time()
mtx = spla.spilu(csc_matrix(B))
end = time.time()
print('Time to decompose B with spilu =', end - start)

计算时间为

Time to decompose B with lu = 4.7765138149261475
Time to decompose B with spilu = 14.165712594985962

1个回答

这种特殊效果很可能来自并行化。

在许多设置中，numpy 将使用多个线程来调用 BLAS/LAPACK 调用。在我的笔记本电脑（Mac OS，本机 Apple python）的默认设置中：

('Time to decompose B with lu =', 9.530492067337036)
('Time to decompose B with spilu =', 20.418880939483643)

活动监视器显示调用期间调用的多个线程，调用期间lu只有一个线程spilu。

在明确指定线程数后（注意，在 Mac 上你也必须这样做）使用“过度杀伤的命令集”：

export MKL_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1
export OMP_NUM_THREADS=1
export VECLIB_MAXIMUM_THREADS=1

香草密集 LU 的时间发生了变化：

('Time to decompose B with lu =', 25.678237199783325)
('Time to decompose B with spilu =', 21.03290104866028)

这个答案也很有用。

一般评论：

scipy.sparse.linalg.spilu对应于稀疏不完全LU 分解，通常用作预条件子。考虑使用专用的稀疏直接求解器，以防它更好地满足您的需求。

其它你可能感兴趣的问题

上一篇如果对其元素的引用次数过多，我应该使用全局堆分配数组还是局部堆栈分配一个？下一篇为提高核心技能而开展的科学计算项目的建议