[SciPy-dev] sparse comments and questions

Peter Skomoroch peter.skomoroch@gmail....
Mon Jul 9 18:33:51 CDT 2007


In Matlab and some C++ libraries, I've used the following "sparse division"
functionality :

When you attempt to divide one sparse matrix by another of the same size,
return the result of elementwise division of the entries in the matrices.
This assumes that the two input matrices have the same sparsity structure.

This is useful when implementing some algorithms like NMF using sparse
matrices.  Right now in scipy, only division of a sparse matrix by a scalar
is supported.

If you look at sparse.py in trunk:

206         def __truediv__(self, other):
207             if isscalarlike(other):
208                 return self * (1./other)
209             else:
210                 raise NotImplementedError, "sparse matrix division not
yet supported"
212         def __div__(self, other):
213             # Always do true division
214             if isscalarlike(other):
215                 return self * (1./other)
216             else:
217                 raise NotImplementedError, "sparse matrix division not
yet supported"

Here is a c implementation of sparse matrix division (mex file which is
called from matlab ... sorry about the formating):



/* spdotdiv.c c = spdotdiv(a,b) Performs matrix element division c=a./b, but
evaluated only at the sparse locations. (a and b must have same sparcity
structure). */

#include "mex.h" #include <string.h> #include <math.h>
#define C (plhs[0]) #define A (prhs[0]) #define B (prhs[1])

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]
) { int m, n, nzmax, nnz; int i; double *apr, *bpr, *cpr;
if (nrhs != 2) mexErrMsgTxt("Two input arguments required."); if
(!mxIsSparse(A) || !mxIsSparse(B)) mexErrMsgTxt("Input arguments must be
m = mxGetM(A); n = mxGetN(A); nzmax = mxGetNzmax(A); nnz = *(mxGetJc(A)+n);
if ((mxGetM(B) != m) || (mxGetN(B) != n) || (mxGetNzmax(B) != nzmax))
mexErrMsgTxt("Input matrices must have same sparcity structure.");
apr = mxGetPr(A); bpr = mxGetPr(B);
if ((C = mxCreateSparse(m,n,nzmax,mxREAL)) == NULL) mexErrMsgTxt("Could not
allocate sparse matrix."); cpr = mxGetPr(C);
memcpy(mxGetIr(C), mxGetIr(A), nnz*sizeof(int)); memcpy(mxGetJc(C),
mxGetJc(A), (n+1)*sizeof(int));
for (i=0; i<nnz; i++) cpr[i] = apr[i]/bpr[i];

Let me know what you think,


On 7/9/07, Nathan Bell <wnbell@gmail.com> wrote:
> On 7/7/07, Peter Skomoroch <peter.skomoroch@gmail.com> wrote:
> > Nathan,
> >
> > Do you have any plans to implement sparse matrix division?  That is
> > something I've found lacking in the sparse matrix support...
> Sorry, I'm not sure what you mean by sparse matrix division.  Can you
> elaborate?
> --
> Nathan Bell wnbell@gmail.com
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev

Peter N. Skomoroch
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-dev/attachments/20070709/1fd34111/attachment.html 

More information about the Scipy-dev mailing list