sparse matrix operations to gain efficiency

[EDITED]

Hello,

I've never really used sparse matrices before but I was wondering if it would be useful in my case to imporve efficiency. My code is as follows

reac=25;spec=20;
kf=rand(1,reac); kb=rand(1,reac);
C=rand(1,spec);
fwd=zeros(reac,spec);bcwd=zeros(reac,spec);

for i=1:reac
    a=randi([1 3],1,5);b=randi([1 spec],1,5);
    fwd(i,b)=a;
    c=randi([1 3],1,5);d=randi([1 spec],1,5);
    bcwd(i,d)=c;
end

i=1:reac;q=kf(i)'.*prod(C.^fwd(i,:),2)-kb(i)'.*prod(C.^bcwd(i,:),2);

sto=fwd-bcwd;

i=1:spec;rate(i)=sum(sto(:,i).*q,1);

I've used an example of reac=25 and spec=20 but in reality I could have reac up to 1500 and spec up to 500. I use the values of rate in differential equations which I solve using ode15s. This is the part of my code that takes the longest, where fwd and bcwd are sparse matrices. Would truning them into sparse matrices using sparse help me? if so, how can I rewrite my code?

An help is appreciated,

Thank you!

[Cedric]

Your code doesn't guarantee that there will be 5 non-zero values per row, as nothing enforces elements of b and d to be unique. Is it the desired behavior?

[AM]

I didn't know how to force the elements of b and d to be unique in my example but since it's maximum 5 non-zero values per row it is not a problem. I also realize now I can simply use the sparse command and the rest of my code does not change.

[Cedric]

It is likely more complex than this.

The gain in performance obtained from using sparse matrices generally involves building vectors I, J, and V of row and column indices and values of non-zero elements, and building a sparse matrix in one shot with a call like

S = sparse( I, J, V, nRows, nCols ) ;

You can easily build I using REPELEM and V using RANDI. But building J is not direct and you must enforce unique values (per row) otherwise the accumlating behavior of SPARSE will produce a side effect that you may not want.

To illustrate this behavior, look at the following:

>> S = sparse( [1 2 1], [1 2 1], [10 20 30], 2, 2 )
S =
   (1,1)       40
   (2,2)       20

Here you can see that as there were two elements (10 and 30) with same indices (1,1), the values were accumulated (summed) and S(1,1)=10+30=40.

PS: if you wanted to enforce uniqueness, you could use RANDPERM, but it won't allow you to generate J in one shot as far as I can tell.

Answer 1

[John D'Errico]

Are those matrices mostly zero? No. It looks like they are 50% zeros. Sparse will not gain anything, because thre is also overhead in working with computations on sparse matrices.

Next, sparse matrices gain when you are doing matrix multiplications, and solving linear systems of eequations, but mainly only then if they are really sparse. 50% zeros does not cut it.

[AM]

Hello,

Sorry, I wrote them like that in that example but in reality they are mostly zero, only having up to 5 non zero elements per row (number of columns can go up to 500).

This is just one part of my code and in other sections I multiply these matrices row by row to other vectors. I just wanted to get an idea of how operations with sparse matrices work but I should have made my example better.

I've edited my code.

[John D'Errico]

Sparse matrices work just like other matrices. So you use them just like you use ANY other normal matrix.

Are you really creating them randomly? If so, then just use the sparse command to make them into sparse matrices. Calling sparse on a full matrix will return a new matrix that is sparse version of that matrix. (The function full does the opposite, as you would expect.)

You could create them directly using the sparse commnd also, but that takes more skill and understanding of how to build a sparse matrix to do it well, and far more than I will teach you.

But you should also recognize that CREATING the sparse matrix as such willl then add some time. You could also find that you won't gain much, if at all. So if some of the time, your matrices are not that sparse, expect it all to take more time when using sparse matrices. Only you can test it to see. But converting to sparse is simple as I said.

[AM]

I realize now there is no difficulty in using sparse matrices here, I will test out my code and see which way is more efficient, thank you for your time!

(I had tried to use sparse matrices but had made a mistake in my code so I thought they were handled differently, sorry)