Эффективно построить матрицу FEM/FVM

Question

Это типичный прецедент для систем уравнений FEM/FVM, поэтому, возможно, это более широкий интерес. С треугольной сеткой порционного

Я хотел бы создать scipy.sparse.csr_matrix. Матричные строки/столбцы представляют значения в узлах сетки. Матрица имеет записи на главной диагонали и где два узла соединены ребром.

Вот MWE, что первый строит node-> edge-> клетки отношения с, а затем строит матрицу:

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. Однако я обнаружил, что это тоже занимает некоторое время. (Не удивительно.)

Еще одна мысль, что у меня было, что я мог бы заменить диагональ V, I, J простым numpy.add.at если была структура csr_matrix -как данных, где хранится главная диагональ в отдельности. Я знаю, что это существует в некоторых других пакетах программного обеспечения, но не может найти его в scipy. Верный?

Возможно, существует разумный способ построения CSR напрямую?

Answer 1

Таким образом, в конце концов, это оказалось, разница между COO-х и КСО sum_duplicates (точно так же, как предположил @hpaulj). Благодаря усилиям всех участников здесь (особенно @ paul-panzer), a PR находится в стадии разработки, чтобы дать tocsr потрясающее ускорение.

SciPy-х tocsr делает lexsort на (I, J), поэтому он помогает организовать индексы таким образом, что (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.

Answer 2

Я бы попытался создать структуру csr напрямую, особенно если вы прибегаете к np.unique, так как это дает вам отсортированные ключи, что составляет половину выполняемой работы.

Я предполагаю, что вы находитесь в точке, где у вас есть i, j сортируется лексикографически и перекрытием v суммируется с помощью np.add.at на дополнительном inverse выходе np.unique.

Тогда 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 

Форсировочная: 5x