[SciPy-dev] Sparse matrix design and broadcasting

[Ed Schofield]

On 07/07/2006, at 7:53 PM, Neilen Marais wrote:

> Hi
>
> I've been looking at the scipy.sparse routines a little, specifically
> the sparse_lil type. Currently fancy indexing seems to be handled by a
> lot of special case code.
>
> I wonder if there is not some way we can use numpy's build in
> broadcasting to make this easier. My thinking is something like this:
>
> 1) lil_matrix.__setitem__() figures out what the equivalent shape  
> of the
> region being assigned to is. Perhaps some numpy routines exist that  
> can
> help here?
>
> e.g.
>
> lil_mat[1,:] -> (1, lil_mat.shape[1])
> lil_mat[[0,1,2], [5,6,7]] -> (3,)
> lil_mat[ix_([0,1,2], [5,6,7])] -> (3,3)
>
> 2) get numpy to broadcast the value being set to that shape
>
> 3) use a single code path that knows how to handle assignation from 2D
> or 1D arrays.
>
> Does this make sense?


Interesting idea. I wouldn't be surprised if we could re-use some  
code from NumPy for this, or perhaps from numarray (whose indexing  
code is written in Python). I went through the numarray code for  
interpreting slices and index arrays before writing the sparse fancy  
indexing code, hoping to re-use whatever I could, but it made quite  
heavy use of strides, and I didn't see how this could be applied to  
sparse matrices. So I more or less reinvented the wheel. But some  
ideas from numarray (in Lib/generic.py) could probably be lifted,  
such as the way the _universalIndexing method handles both getting  
and setting.

Maybe we could indeed use functionality from NumPy in interpreting  
indices, but I don't know whether this is possible or how we'd do it.  
I'd appreciate any help.

-- Ed