Primitive operations on sparse Matrix done on full copy??

Primitive operations on sparse Matrix done on full copy??

I have a number of sparse matrices, that I want to sum up. These are
composed in a for loop, and summing them up leeds to memory allocation
and reordering in every step of the loop. Very ugly and slow. So a
fellow came up with the idea of putting these matrices in columns of
large sparse matrix and do a

sum(M,2).

This leeds to a out of memory error. A

sum(M)

works fine. So I assume, Matlab does the summing along rows on
full(M),
wich is easier to implement, but not of that much use, as in the above
stated example, one dimension is some 10 million.
Yes, you read right, it is very sparse.

Transposing M seems to be done on full(M) as well.

I tried to do the sum on every row at a time, but even on row is too
large (otherwise the original matrices wouldn't have to be sparse ;-))

I ended up with doing the sum by hand, but a for loop over every
nonezero of the final matrix is REALLY slow.

Any good workaround?

(to Mathworks: if it isn't yet, would you copy these 'features' from my
wishlist to yours? Thanx)

Axel

Primitive operations on sparse Matrix done on full copy??

Quote:> I have a number of sparse matrices, that I want to sum up. These are
> composed in a for loop, and summing them up leeds to memory allocation
> and reordering in every step of the loop. Very ugly and slow. So a
> fellow came up with the idea of putting these matrices in columns of
> large sparse matrix and do a

> sum(M,2).

> This leeds to a out of memory error. A

> sum(M)

> works fine. So I assume, Matlab does the summing along rows on
> full(M),
> wich is easier to implement, but not of that much use, as in the above
> stated example, one dimension is some 10 million.
> Yes, you read right, it is very sparse.

> Transposing M seems to be done on full(M) as well.

> I tried to do the sum on every row at a time, but even on row is too
> large (otherwise the original matrices wouldn't have to be sparse ;-))

> I ended up with doing the sum by hand, but a for loop over every
> nonezero of the final matrix is REALLY slow.

> Any good workaround?

Hi.

If A is an m-by-n sparse matrix with nz = nnz(A) nonzero elements,
the computational complexity of these kinds of primitive operations
is proportional to nz and to n, but is pretty much independent of m.
We never do a sparse matrix operation on A by forming full(A) unless
the final result requires that much memory.  So the memory requirements
and execution time are never proportional to m*n unless the result
requires that.  In other words, if your sparse matrix is highly
rectangular, orient it so that it is tall and skinny, not short and fat.

I'm afraid your description doesn't provide enough detail for me to
offer any more specific advice.  If you want to contact me, or

code, we'll try to help.

-- Cleve Moler

Primitive operations on sparse Matrix done on full copy??

Quote:>I have a number of sparse matrices, that I want to sum up. These are
>composed in a for loop, and summing them up leeds to memory allocation
>and reordering in every step of the loop. Very ugly and slow. So a
>fellow came up with the idea of putting these matrices in columns of
>large sparse matrix and do a

>sum(M,2).

>This leeds to a out of memory error. A

>sum(M)

>works fine. So I assume, Matlab does the summing along rows on
>full(M),
>wich is easier to implement, but not of that much use, as in the above
>stated example, one dimension is some 10 million.

I doubt that Matlab is that stupid (at least Matlab5), but it's true
that sum(M,2) is very inefficient, but not _that_ inefficient.

I would guess that the reason is that Matlab uses shiftdim to rearrange
the dimensions, and for some reasons shiftdim is not fast.

A simple solution would be to use [I,J,S]=find(M);sparse(I,1,S); which
generates the right answer, but it is also too slow.

A better solution is to write a c-mex function for the specific task,
as follows: (save as sum2.c and compile with mex -O sum2.c).
----

/* Computes sum(M,2)  fast for a sparse matrix */
#include <math.h>
#include "mex.h"

void mexFunction(int nlhs,mxArray*plhs[],int nrhs,const mxArray*prhs[])
{
int n,nnz,i;
const double *s;
const int *index;
double *out;
if (nrhs!=1) mexErrMsgTxt("sum2.c: Exactly 1 argument required");
if (!mxIsSparse(prhs)) mexErrMsgTxt("sum2.c: Only for sparse matrices");
n=mxGetM(prhs);
nnz=mxGetJc(prhs)[mxGetN(prhs)];
plhs=mxCreateDoubleMatrix(n,1,mxREAL);
out=mxGetPr(plhs);
index=mxGetIr(prhs);
s=mxGetPr(prhs);
for(i=0;i<nnz;i++,index++,s++) {
int in=*index;
out[in]+=*s;
};

Quote:}

----
Testing it on 1e5*1e5 random sparse matrices with approximately 1e6 nnz
elements, gave this (all approximate and computed on 64M ultra-1):

sprandn(1e5,1e5,1e-4)    330 seconds
sum(M,2)                 70 seconds
M'                       57 seconds
[I,J,S]=find(M)         101 seconds
sum(M)                   1-14 seconds depending on caching
sum2(M)                0.5-24 seconds depending on caching

So the new C-function is as fast as traditional sum, but does what you want.

I don't know why find and transpose are so slow, because
it seems as though find only has to do one pass through the matrix and
then copy elements to two arrays and easily create a third.

Transpose could be written as at most two passes through the matrix,
which would take a few seconds.
--
// Homepage  http://www.dna.lth.se/home/Hans_Olsson/

Primitive operations on sparse Matrix done on full copy??

This was the thing I was hoping for. So thanx to Hans Olsson.

But still I'm not satisfied, sorry. This mex-file returns a full vector.
But in my case still that vector is sparse, so did anybody (Hans??) do
this with sparse output??
I did a quick check and it get's a bit more complicated, so if anybody
did this, please say so.

Otherwise I volunteer to do my learning in mex-files on this and will
post the result (hopefully :-))

Axel

Primitive operations on sparse Matrix done on full copy??

>>I have a number of sparse matrices, that I want to sum up. These are
>>composed in a for loop, and summing them up leeds to memory allocation
>>and reordering in every step of the loop. Very ugly and slow. So a
>>fellow came up with the idea of putting these matrices in columns of
>>large sparse matrix and do a

>>sum(M,2).

>>This leeds to a out of memory error. A

>>sum(M)

>>works fine. So I assume, Matlab does the summing along rows on
>>full(M),
>>wich is easier to implement, but not of that much use, as in the above
>>stated example, one dimension is some 10 million.

[Deleted]

Quote:>A simple solution would be to use [I,J,S]=find(M);sparse(I,1,S); which
>generates the right answer, but it is also too slow.

After posting this answer, I found out that I'd forgotten the most
obvious solution:

The simplest solution is of course X*ones(size(X,2),1) which actually
is several times faster than sum(X,2), and almost as fast as the c-code
I gave.

I think this clearly shows that matlab's routines for sparse matrices
are not optimal (summing the elements should not be slower than multiplying
by one and then summing them).

Quote:>Testing it on 1e5*1e5 random sparse matrices with approximately 1e6 nnz
>elements, gave this (all approximate and computed on 64M ultra-1):

>sprandn(1e5,1e5,1e-4)    60-330 seconds depending on caching
>sum(M,2)                 70 seconds
>M'                       57 seconds
>[I,J,S]=find(M)         101 seconds
>sum(M)                   1-14 seconds depending on caching
>sum2(M)                0.5-24 seconds depending on caching

X*ones(1e5,1)            0.6-30 seconds depending on caching

Quote:>I don't know why find and transpose are so slow, because
>it seems as though find only has to do one pass through the matrix and
>then copy elements to two arrays and easily create a third.

>Transpose could be written as at most two passes through the matrix,
>which would take a few seconds.

Experimenting showed that the reason for the reason for the slow find and
transpose seems to be the bad use of caches etc. Basically the data hardly
fits into memory and then performancy drops dramatically. I'm not sure
how performance related to memory can be fixed, and my experience with
Matlab is that the algorithms are not written to avoid cache problems.

--
// Homepage  http://www.dna.lth.se/home/Hans_Olsson/

Primitive operations on sparse Matrix done on full copy??

Hi there,
I volunteered and succeeded.
I have rewritten the mex-file given by Hans to output a sparse vector.

This was generally not that hard, but there is one tricky thing I like
to mention, as it is subject to some optimization: how do I get the
index to pr for a given row. I do this by guessing a starting point and
then just walk the list, which is quite good for equally distibuted
matrices (sum-vector to be precise). It is easy done and fairly good.

I'd be glad if Hans could include this listing in his benchmarks.

I ran a few tests on small matrices and the long benchmark is running
right now, not completed yet, but I'll leave office soon.

Anyway, here it goes, if there are hints for optimization, I'm glad to
(remember, its my first (half) mex-file :-)

Enjoy...

Axel
---------------------------

/* Computes sum(M,2)  fast for a sparse matrix */
/* Modified by A. Hecht for producing sparse output */

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "mex.h"

static int cmpval(const void* i1,const void* j2)
{
int *i,*j;
i=(int*)i1;
j=(int*)j2;
if (*i>*j)
return (1);
if (*i<*j)
return (-1);
return(0);

Quote:}

void mexFunction(int nlhs,mxArray*plhs[],int nrhs,const mxArray*prhs[])
{
int n,nnz,i,in,nunik;
const double *s;
const int *index;
int *field, *unik, old;
double *out;
if (nrhs!=1) mexErrMsgTxt("sum2.c: Exactly 1 argument required");
if (!mxIsSparse(prhs)) mexErrMsgTxt("sum2.c: Only for sparse
matrices");
n=mxGetM(prhs);
nnz=mxGetJc(prhs)[mxGetN(prhs)];
index=mxGetIr(prhs);
s=mxGetPr(prhs);
field = mxCalloc(nnz,sizeof(int));
unik = mxCalloc(nnz,sizeof(int));
memmove(field,index,nnz*sizeof(int));
qsort(field,nnz,sizeof(int),cmpval);
old=field;
nunik=1;
for (i=1;i<nnz;i++)
{
if (old!=field[i])
{
old=field[i];
unik[nunik++]=field[i];
}
};
plhs=mxCreateSparse(n,1,nunik,mxREAL);
out=mxGetPr(plhs);
memmove(mxGetIr(plhs),unik,nunik*sizeof(int));
for(i=0;i<nnz;i++)
{
in=floor(nunik*index[i]/nnz);
while ( unik[in] != index[i])
{
if (unik[in]>index[i])
{
in--;
}
else
{
in++;
}
};
out[in]+=s[i];
};
mxGetJc(plhs)=nunik;
Quote:}

Primitive operations on sparse Matrix done on full copy??

Sorry folks, this is the errata:
I forgot the initial entry in ir of the output matrix (my test matrices
had a result in the first row, so that error didn't show :-( )

And there was an error, which I had corrected, but it reappeared on some
strange path (which may be my brain or the undo in emacs and propably a
combination of the two).
This is the guessed starting value in the sum-loop, which is now
corrected and multiplied by 1.0 to avoid overflows (yielding strange
errors, as one might expect)

So here comes the working version:

-------------
/* Computes sum(M,2)  fast for a sparse matrix */
/* Modified by A. Hecht for producing sparse output */

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "mex.h"

static int cmpval(const void* i1,const void* j2)
{
int *i,*j;
i=(int*)i1;
j=(int*)j2;
if (*i>*j)
return (1);
if (*i<*j)
return (-1);
return(0);

Quote:}

void mexFunction(int nlhs,mxArray*plhs[],int nrhs,const mxArray*prhs[])
{
int n,nnz,i,in,nunik;
const double *s;
const int *index;
int *field, *unik, old;
double *out;
if (nrhs!=1) mexErrMsgTxt("sum2.c: Exactly 1 argument required");
if (!mxIsSparse(prhs)) mexErrMsgTxt("sum2.c: Only for sparse
matrices");
n=mxGetM(prhs);
nnz=mxGetJc(prhs)[mxGetN(prhs)];
index=mxGetIr(prhs);
s=mxGetPr(prhs);
field = mxCalloc(nnz,sizeof(int));
unik = mxCalloc(nnz,sizeof(int));
memmove(field,index,nnz*sizeof(int));
qsort(field,nnz,sizeof(int),cmpval);
old=field;
unik=old;  /*errata 1*/
nunik=1;
for (i=1;i<nnz;i++)
{
if (old!=field[i])
{
old=field[i];
unik[nunik++]=field[i];
}
};
plhs=mxCreateSparse(n,1,nunik,mxREAL);
out=mxGetPr(plhs);
memmove(mxGetIr(plhs),unik,nunik*sizeof(int));
for(i=0;i<nnz;i++)
{
in=floor(1.0*nunik*index[i]/n);   /*errata 2*/
while ( unik[in] != index[i])
{
if (unik[in]>index[i])
{
in--;
}
else
{
in++;
}
};
out[in]+=s[i];
};
mxGetJc(plhs)=nunik;
Quote:}

Hi,

I found that there is a difference of sparse matrix operations between
version 4.x and 5.x.  If X is a sparse matrix, in version 4.x I can do
something like
Y(1:20,:) = X(5:24,:)
.  (This is just a simplified example of a more complicated task.)
However, when I try to do the same with version 5.x I got the message
"???  In an assignment  A(:,matrix) = B, the number of rows in A   and B
must be the same."
It seems I can only get around by first convert the sparse matrix to
full matrix, X = full(X), and then do the operation.  It seems to me
that this approach violates the purpose of having sparse matrix in the
first place, and is not practical to me since the matrix X is very
large.  Is there anyway I can preserve X as a sparse matrix and still
accomplish the operation?  Thanks in advance.

H.J. Wang