CODEWITHSUNDEEP
pythonnumpyscipysparse-matrix
Lin
Lin

Reputation: 1203

How to avoid sparse to dense matrix convertions

I have a simple code that I work with sparse matrices in Numpy/Scipy the code is shown below:

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()))))

I would like to remove the sparse to dense convertions, but dont know how to replace the toeplitz function and the fliplr withtout the convertion. Right now I have this:

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

Of course I can work with non sparse matrices and convert just in the end, but I'm looking to work always with sparse ones. Any better ideas?

Upvotes: 0

Views: 275

Answers (1)

Paul Panzer
Paul Panzer

Reputation: 53029

Here is how to do it using scipy.sparse.diags. diags is similar to spdiags but a bit more convenient.

>>> 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))

This example has row vectors, for column vectors one can cast to csc and then proceed in the same way.

>>> # 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 can be done by indexing. Small gotcha: not all sparse matrix classes do currently support indexing, you may have to cast.

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

Upvotes: 1

Related Questions