高效构建 FEM/FVM 矩阵答案

【问题标题】：Efficiently construct FEM/FVM matrix高效构建 FEM/FVM 矩阵
【发布时间】：2017-02-22 15:56:20
【问题描述】：

这是 FEM/FVM 方程系统的典型用例，因此可能具有更广泛的兴趣。从三角形网格 à la

我想创建一个scipy.sparse.csr_matrix。矩阵行/列表示网格节点处的值。该矩阵在主对角线上以及两个节点通过边连接的任何位置都有条目。

这是一个 MWE，它首先建立一个节点->边->细胞关系，然后建立矩阵：

import numpy
import meshzoo
from scipy import sparse

nx = 1600
ny = 1000
verts, cells = meshzoo.rectangle(0.0, 1.61, 0.0, 1.0, nx, ny)

n = len(verts)

nds = cells.T
nodes_edge_cells = numpy.stack([nds[[1, 2]], nds[[2, 0]],nds[[0, 1]]], axis=1)

# assign values to each edge (per cell)
alpha = numpy.random.rand(3, len(cells))
vals = numpy.array([
    [alpha**2, -alpha],
    [-alpha, alpha**2],
    ])

# Build I, J, V entries for COO matrix
I = []
J = []
V = []
#
V.append(vals[0][0])
V.append(vals[0][1])
V.append(vals[1][0])
V.append(vals[1][1])
#
I.append(nodes_edge_cells[0])
I.append(nodes_edge_cells[0])
I.append(nodes_edge_cells[1])
I.append(nodes_edge_cells[1])
#
J.append(nodes_edge_cells[0])
J.append(nodes_edge_cells[1])
J.append(nodes_edge_cells[0])
J.append(nodes_edge_cells[1])
# Create suitable data for coo_matrix
I = numpy.concatenate(I).flat
J = numpy.concatenate(J).flat
V = numpy.concatenate(V).flat

matrix = sparse.coo_matrix((V, (I, J)), shape=(n, n))
matrix = matrix.tocsr()

与

python -m cProfile -o profile.prof main.py
snakeviz profile.prof

可以创建和查看上述个人资料：

tocsr() 方法在这里占据了运行时的最大份额，但在构建 alpha 更复杂时也是如此。因此，我正在寻找加快速度的方法。

我已经找到了什么：

由于数据的结构，矩阵对角线上的值可以提前求和，即

V.append(vals[0, 0, 0] + vals[1, 1, 2])
I.append(nodes_edge_cells[0, 0])  # == nodes_edge_cells[1, 2]
J.append(nodes_edge_cells[0, 0])  # == nodes_edge_cells[1, 2]

这使得I、J、V 更短，从而加快了tocsr。

现在，边缘是“每个单元格”的。我可以使用numpy.unique 识别彼此相等的边缘，有效节省了大约一半的I、J、V。但是，我发现这也需要一些时间。（不足为奇。）

我的另一个想法是，如果有一个类似于csr_matrix 的数据结构，其中主对角线为分开存放。我知道这存在于其他一些软件包中，但在 scipy 中找不到。对吗？

也许有一种明智的方式来直接构建 CSR？

【问题讨论】：

我想看看不同步骤的一些时间安排。 tocsr 使用编译后的代码。我认为事先对coo 输入的任何按摩都需要同样长的时间。
您确定tocsr 需要很长时间吗？我对一个 10k x 10k 矩阵做了非常相似的事情，其中 I、J、V 的长度达到了数百万，而且不需要那么长时间。也许 5-10 秒。
感谢 cmets。我已在原始帖子中添加了个人资料。
有趣。我认为sum_duplicates 是编译后的csr 构造的一部分，但您的配置文件显示它是coo.py 中的Python 方法。它使用lexsort、np.nonzero 和reduceat。所以做你自己的sum_duplicates 是有希望的。如果你的 coo has_canonical_format tocsr 应该更快。
csr 有自己的sum_duplicates 方法。它使用编译后的sparse._sparsetools.csr_sum_duplicates。

标签： python numpy matrix scipy scientific-computing

【解决方案1】：

我会尝试直接创建 csr 结构，特别是如果您使用 np.unique，因为这会为您提供排序的键，这是完成了一半的工作。

我假设您处于 i, j 按字典顺序排序并在 np.unique 的可选 inverse 输出上使用 np.add.at 对 v 求和的位置。

那么v和j已经是csr格式了。剩下要做的就是创建indptr，您只需通过np.searchsorted(i, np.arange(M+1)) 获得，其中M 是列长度。您可以将这些直接传递给sparse.csr_matrix 构造函数。

好的，让代码说话：

import numpy as np
from scipy import sparse
from timeit import timeit


def tocsr(I, J, E, N):
    n = len(I)
    K = np.empty((n,), dtype=np.int64)
    K.view(np.int32).reshape(n, 2).T[...] = J, I  
    S = np.argsort(K)
    KS = K[S]
    steps = np.flatnonzero(np.r_[1, np.diff(KS)])
    ED = np.add.reduceat(E[S], steps)
    JD, ID = KS[steps].view(np.int32).reshape(-1, 2).T
    ID = np.searchsorted(ID, np.arange(N+1))
    return sparse.csr_matrix((ED, np.array(JD, dtype=int), ID), (N, N))

def viacoo(I, J, E, N):
    return sparse.coo_matrix((E, (I, J)), (N, N)).tocsr()


#testing and timing

# correctness
N = 1000
A = np.random.random((N, N)) < 0.001
I, J = np.where(A)

E = np.random.random((2, len(I)))
D = np.zeros((2,) + A.shape)
D[:, I, J] = E
D2 = tocsr(np.r_[I, I], np.r_[J, J], E.ravel(), N).A

print('correct:', np.allclose(D.sum(axis=0), D2))

# speed
N = 100000
K = 10

I, J = np.random.randint(0, N, (2, K*N))
E = np.random.random((2 * len(I),))
I, J, E = np.r_[I, I, J, J], np.r_[J, J, I, I], np.r_[E, E]

print('N:', N, ' --  nnz (with duplicates):', len(E))
print('direct: ', timeit('f(a,b,c,d)', number=10, globals={'f': tocsr, 'a': I, 'b': J, 'c': E, 'd': N}), 'secs for 10 iterations')
print('via coo:', timeit('f(a,b,c,d)', number=10, globals={'f': viacoo, 'a': I, 'b': J, 'c': E, 'd': N}), 'secs for 10 iterations')

打印：

correct: True
N: 100000  --  nnz (with duplicates): 4000000
direct:  7.702431229001377 secs for 10 iterations
via coo: 41.813509466010146 secs for 10 iterations

加速：5 倍

【讨论】：

我实际上想避免使用unique，因为它不是 O(n) 操作（或者我认为）。当然，如果to_csr 还是这样做了，我还不如自己做。所有好的提示。
我似乎记得 unique 使用 argsort，所以它是 O(n log n)。另一方面，CSR 对索引进行了排序，因此除非您的数据以任何形式预先排序（是吗？）或者您有一些图论约束（例如，每个节点的边数固定或有界），否则您无法避免 O(n日志 n)。我不认为您的数据无论如何都按节点或类似的东西分组？网格是否像看起来一样规则？因为可以利用这种结构进行更快的排序。
看起来加快速度的关键是以最小化lexsort 时间的方式生成元素。如果重复项彼此相邻，add.reduceat 可以快速求和。
@NicoSchlömer 我已经实现了一个直接的 tocsr，它似乎比通过 coo 快很多。很可能大部分的加速是由于我避免lexsort 支持argsort。你能用你的实际数据检查一下吗？我很想知道它在现实世界中的效果如何。
@hpaulj 似乎 lexsort 使用 mergesort 而 argsort 默认使用 quicksort。当我通过kind=mergesort 时，argsort 产生与 lexsort 相同的输出。并且运行速度慢了两倍多。但仍然是 lexsort 的两倍多。这些差异确实令人困惑。

【解决方案2】：

所以，最终结果证明这是首席运营官和企业社会责任 sum_duplicates 之间的区别（就像 @hpaulj 怀疑的那样）。感谢所有参与其中的人（尤其是@paul-panzer）的努力，a PR 正在为tocsr 提供巨大的加速。

SciPy 的tocsr 在(I, J) 上执行lexsort，因此它有助于以(I, J) 已经得到公平排序的方式组织索引。

对于上例中的nx=4、ny=2，I 和J 是

[1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7]
[1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7 1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7]

先对cells的每一行排序，再按第一列的行排序

cells = numpy.sort(cells, axis=1)
cells = cells[cells[:, 0].argsort()]

生产

[1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6]
[1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6 1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6]

对于原始帖子中的数字，排序将运行时间从大约 40 秒缩短到 8 秒。

如果一开始就对节点进行更适当的编号，也许可以实现更好的排序。我在想Cuthill-McKee 和friends。

【讨论】：