Как избежать редких преобразований в плотные матрицы - PullRequest
Новая
       11

Как избежать редких преобразований в плотные матрицы

0 голосов
/ 07 октября 2018

У меня есть простой код, который я работаю с разреженными матрицами в Numpy / Scipy, код показан ниже: 

import numpy as np
import scipy as sp
from scipy.sparse import csr_matrix as sparse
from scipy.sparse import vstack
from scipy.linalg import toeplitz

N = 100
d4 = np.array([-22/24, 17/24, 9/24 ,-5/24, 1/24])
s4 = np.array([1/24,-9/8,9/8,-1/24])
n = len(s4)
r4 = sparse((s4, (np.zeros(n), np.arange(n))), shape=[1, N+1])
c4 = sparse(([s4[0]], ([0], [0])), shape=[N-2, 1])
lnd = len(d4)
rd4 = sparse((d4, (np.zeros(lnd), np.arange(lnd))), shape=[1, N+1])
D = sparse(np.concatenate((rd4.todense(), toeplitz(c4.todense(),r4.todense()), np.fliplr(rd4.todense()))))

Я хотел бы удалить преобразования sparse в dense, но незнать, как заменить функцию toeplitz и fliplr без преобразования.Прямо сейчас у меня есть это: 

D = vstack([rd4, sparse(toeplitz(c4.todense(),r4.todense())), sparse(np.fliplr(rd4.todense()))])

Конечно, я могу работать с не разреженными матрицами и конвертировать только в конце, но я всегда стараюсь работать с разреженными.Есть идеи получше?

1 Ответ

0 голосов
/ 07 октября 2018

Вот как это сделать, используя scipy.sparse.diags.diags похож на spdiags, но немного удобнее. 

>>> import numpy as np
>>> from scipy import sparse, linalg
>>> 
>>> a = sparse.csr_matrix(np.random.randint(-50, 10, (1, 10)).clip(0, None))
>>> b = sparse.csr_matrix(np.random.randint(-50, 10, (1, 10)).clip(0, None))

В этом примере есть векторы строк, для векторов столбцов можно привести к csc, а затем продолжить таким же образом. 

>>> # dense method for reference
>>> d_toepl = linalg.toeplitz(a.A, b.A)
>>> 
>>> idx = b.indices[0] == 0 # make sure first element of b is ignored
>>> vals, offs = np.r_[a.data, b.data[idx:]], np.r_[-a.indices, b.indices[idx:]]
>>> N = max(a.shape[1], b.shape[1])
>>> dtype = (a[0, ...] + b[0, ...]).dtype
>>> 
>>> s_toepl = sparse.diags(vals, offs, (N, N), dtype=dtype)
>>> 
>>> np.all(d_toepl == s_toepl)
True

fliplr можно сделать с помощью индексации.Небольшое замечание: не все классы разреженных матриц в настоящее время поддерживают индексирование, возможно, вам придется приводить. 

>>> np.all(np.fliplr(d_toepl) == s_toepl.tocsr()[:, ::-1])
True
...