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?

Solution

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