Преобразование из разреженного в плотное в разреженное снова уменьшает плотность после построения разреженной матрицы - PullRequest
Новая
       15

Преобразование из разреженного в плотное в разреженное снова уменьшает плотность после построения разреженной матрицы

2 голосов
/ 14 марта 2019

Я использую scipy для генерации разреженной матрицы конечных разностей, сначала строю ее из блочных матриц, а затем редактирую диагональ для учета граничных условий. Полученная разреженная матрица имеет тип BSR. Я обнаружил, что если я преобразую матрицу в плотную матрицу, а затем обратно в разреженную матрицу с помощью функции scipy.sparse.BSR_matrix, у меня останется более разреженная матрица, чем раньше. Вот код, который я использую для генерации матрицы: 

size = (4,4)

xDiff = np.zeros((size[0]+1,size[0]))
ix,jx = np.indices(xDiff.shape)
xDiff[ix==jx] = 1
xDiff[ix==jx+1] = -1

yDiff = np.zeros((size[1]+1,size[1]))
iy,jy = np.indices(yDiff.shape)
yDiff[iy==jy] = 1
yDiff[iy==jy+1] = -1

Ax = sp.sparse.dia_matrix(-np.matmul(np.transpose(xDiff),xDiff))
Ay = sp.sparse.dia_matrix(-np.matmul(np.transpose(yDiff),yDiff))

lap = sp.sparse.kron(sp.sparse.eye(size[1]),Ax) + sp.sparse.kron(Ay,sp.sparse.eye(size[0]))

#set up boundary conditions
BC_diag = np.array([2]+[1]*(size[0]-2)+[2]+([1]+[0]*(size[0]-2)+[1])*(size[1]-2)+[2]+[1]*(size[0]-2)+[2])

lap += sp.sparse.diags(BC_diag)

Если я проверяю разреженность этой матрицы, я вижу следующее: 

lap
<16x16 sparse matrix of type '<class 'numpy.float64'>'
with 160 stored elements (blocksize = 4x4) in Block Sparse Row format>

Однако, если я преобразую его в плотную матрицу, а затем вернусь в тот же разреженный формат, я увижу гораздо более разреженную матрицу: 

sp.sparse.bsr_matrix(lap.todense())
<16x16 sparse matrix of type '<class 'numpy.float64'>'
with 64 stored elements (blocksize = 1x1) in Block Sparse Row format>

Я подозреваю, что причина этого заключается в том, что я построил матрицу, используя функцию sparse.kron, но мой вопрос заключается в том, есть ли способ получить меньшую разреженную матрицу, не преобразуя сначала в плотную, например, если я закончу желая смоделировать очень большой домен.

Ответы [ 2 ]

1 голос
/ 14 марта 2019

BSR сохраняет данные в плотных блоках: 

In [167]: lap.data.shape                                                        
Out[167]: (10, 4, 4)

В этом случае эти блоки имеют довольно много нулей. 

In [168]: lap1 = lap.tocsr() 
In [170]: lap1                                                                  
Out[170]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 160 stored elements in Compressed Sparse Row format>
In [171]: lap1.data                                                             
Out[171]: 
array([-2.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,  1., -3.,  1.,  0.,  0.,
        1.,  0.,  0.,  0.,  1., -3.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,
        1., -2.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0., -3.,  1.,  0.,
        0.,  1.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  1., -4.,  1.,  0., 
        ...
        0.,  0.,  1., -2.])

Очистка на месте: 

In [172]: lap1.eliminate_zeros()                                                
In [173]: lap1                                                                  
Out[173]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 64 stored elements in Compressed Sparse Row format>

Если указать формат csr при использовании kron: 

In [181]: lap2 = sparse.kron(np.eye(size[1]),Ax,format='csr') + sparse.kron(Ay,n
     ...: p.eye(size[0]), format='csr')                                         
In [182]: lap2                                                                  
Out[182]: 
<16x16 sparse matrix of type '<class 'numpy.float64'>'
    with 64 stored elements in Compressed Sparse Row format>
0 голосов
/ 14 марта 2019

[Мне сообщили, что мой ответ неверен.Причина, если я понимаю, в том, что Scipy не использует Lapack для создания матриц, а использует для этой цели свой собственный код.Интересно.Информация, хотя и неожиданная, имеет авторитетный характер.Я буду откладывать на это!

[Я оставлю ответ опубликованным для справки, но больше не буду утверждать, что ответ был правильным.]

Вообще говоря, когда речь идет о сложных структурах данных, таких какразреженные матрицы, у вас есть два случая:

  1. конструктор заранее знает все содержимое структуры;или
  2. структура спроектирована для постепенного построения, чтобы полное содержание структуры стало известно только после ее завершения.

Классическим случаем сложной структуры данных являетсяслучай двоичного дерева.Вы можете сделать бинарное дерево более эффективным, скопировав его после его завершения.В противном случае стандартная красно-черная реализация дерева оставляет некоторые пути поиска в два раза длиннее других - что обычно нормально, но не оптимально.

Теперь вы, вероятно, знаете все это, но я упоминаюэто по причине.Сципи зависит от Лапака.Lapack предлагает несколько разных схем хранения.Двумя из них являются

  • общие разреженные и
  • полосатые

схемы.Похоже, что Scipy начинает с того, что хранит вашу матрицу как разреженную, где индексы каждого ненулевого элемента хранятся в явном виде;но при копировании Сципи замечает, что полосатое представление является более подходящим - ведь ваша матрица, в конце концов, полосатая.

...