在Python中如何快速组装稀疏矩阵

在采用Python编写数值计算方法的求解器时，往往需要自己组装求解的刚度矩阵，当矩阵带宽比较大时，采用稠密矩阵（尽管直接组装速度会快些）常常会遇到内存受限的问题。此时，我们会考虑采用稀疏矩阵来存储，这样可以省去为零位置的值，大大节省内存或磁盘的占用。
在Python中组装刚度矩阵时，这里采用NumPy以及SciPy两个科学计算依赖库。

NumPy

NumPy是Python中科学计算的基本软件包。它是一个Python库，提供多维数组对象，各种派生对象（例如蒙版数组和矩阵）以及各种例程，用于对数组进行快速操作，包括数学，逻辑，形状处理，排序，选择，I/O ，离散傅立叶变换，基本线性代数，基本统计运算，随机模拟等等。

在安装时需注意其Fortran或C++依赖库laspack、bla等是否存在。推荐使用conda进行安装，其默认采用Intel数学核心库进行编译。MKL较OpenBLAS更快以及更稳定。

SciPy

数值算法和特定领域的工具箱，包括信号处理，优化，统计等。
这里仅使用到Sparse linear algebra模块

在组装过程中，采用SciPy中的Coo_Matrix构造稀疏矩阵，只需构造row，col以及data。在构造过程中，存在三个难点，在此记录。

提取ndarray中的非零元素

比较Python中常用的几种方式的提取效率

import timeit
import numpy as np
import cv2
def check_speed(a,positions):
    print('\t Number of nonzero elements:',np.sum(a>0))
    print('\t Using cv2:', timeit.timeit(lambda:cv2.findNonZero((a)).squeeze()[:,0],number=100)/100.0, 's')
    print('\t Using numpy nonzero:', timeit.timeit(lambda:np.nonzero(a)[1],number=100)/100.0, 's')
    print('\t Using readily made matrix:', timeit.timeit(lambda:positions[(np.abs(a)>1e-16)[0]],number=100)/100.0 ,'s')
    print('\t Using readily made matrix:', timeit.timeit(lambda:np.arange(a.size),number=100)/100.0 ,'s')

output

Number of nonzero elements: 499413
Using cv2: 0.02275758675998077 s
Using numpy nonzero: 0.007664396830368787 s
Using readily made matrix: 0.0030909083294682204 s

可以看出使用readily made matrix的方式速度最快，因为这里处理的是一维数组，对于多维数组而言，cv2为最优，numpy nonzero次之。

组装row，col，data

遍历刚度矩阵中的每一行，将数据填入row，col，data数组中，我们一般采用NumPy中的concatenate方法进行连接，值得注意的是，不应该一边遍历一边连接，最有效的方法应该是将每次遍历的结果Append至列表中，最后采用concatenate一次性连接。

如何让节点的遍历更快速

一般而言，遍历每一个点的操作应该是相互独立的，可对这部分进行并行编程，充分利用CPU的多核性能进行计算，Python中可考虑采用multiprocessing模块进行处理。当然，这当中，代码的设计、中间量的内存管理等至关重要。

参考