如何在Python/Scipy中有效地组装大型稀疏矩阵

Question

3

我正在使用Scipy进行FEM项目的工作。现在我的问题是，稀疏矩阵的组装速度太慢了。我计算每个元素在密集小矩阵中的贡献（每个元素一个小矩阵）。为了组装全局矩阵，我循环遍历所有小密集矩阵，并按以下方式设置矩阵条目：

[i,j] = someList[k][l]
Mglobal[i,j] = Mglobal[i,j] + Mlocal[k,l]

Mglobal是适当大小的lil_matrix，someList映射索引变量。

当然，这种方法相当慢，并且消耗了大部分矩阵组装时间。有没有更好的方法可以从许多小而密集的矩阵中组装大稀疏矩阵？我尝试使用scipy.weave，但似乎它不适用于稀疏矩阵。

- cp3028

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- David Wells · Accepted Answer

我在scipy邮件列表上发布了我的回复; 访问stack overflow更容易，因此我将稍微改进的版本也发布在这里。

技巧是使用IJV存储格式。这是三个数组的三元组，第一个包含行索引，第二个包含列索引，第三个包含该位置矩阵的值。这是构建有限元矩阵（或任何稀疏矩阵，在我看来）的最佳方法，因为访问此格式非常快速（只需填充数组）。

在scipy中，这称为coo_matrix; 该类将三个数组作为参数。它实际上只有用于转换为另一种格式（CSR os CSC），以进行快速线性代数运算。

对于有限元素，您可以通过以下方式估计三个数组的大小。

size = number_of_elements * number_of_basis_functions**2

如果您有2D二次方程，例如，您将执行number_of_elements * 36。这种方法很方便，因为如果您具有本地矩阵，则绝对具有全局数字和条目值：正是构建三个IJV数组所需的内容。 Scipy足够聪明，可以抛出零条目，因此高估没有问题。