[SciPy-dev] sparse comments and questions

[Peter Skomoroch]

Nathan,

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"
211
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):

source:

http://journalclub.mit.edu/jclub/message?com_id=2;publication_id=21;message_id=58;session_id=2E87004B582D814FFB7D7DC3E64C3789;seq_no=54958

/* 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
sparse.");
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,

-Pete

On 7/9/07, Nathan Bell <wnbell at gmail.com> wrote:
>
> On 7/7/07, Peter Skomoroch <peter.skomoroch at 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 at gmail.com
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>



-- 
Peter N. Skomoroch
peter.skomoroch at gmail.com
http://www.datawrangling.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20070709/1fd34111/attachment.html>