 good chemistry is complicated,and a little bit messy -LW PerlMonks

### Sparse Matrix Multiplication Problem

by neversaint (Deacon)
 on Jan 27, 2009 at 11:24 UTC Need Help??

neversaint has asked for the wisdom of the Perl Monks concerning the following question:

Dear Masters,
Given this sparse matrix representation:
```my @realval = (
0.994704478,
0.989459074,
0.994717023,
1.000000000,
1.000000000,
0.002647761,
0.005282977,
0.000882587,
0.005270463,
0.000882587,
0.002635231,
0.000882587,
0.002635231,);

my @row_index = qw( 1 2 3 4 5 1 3 1 2 1 2 1 2);
my @col_index =  qw(1 2 3 4 5 5 3 2 1 3 3 4 4);
my @mat_dim  = (5,5); # M always equal to N
I have no problem generating the actual matrix like this:
``` my @actual_mat = (

[ 0.994704478, 0.000882587, 0.000882587, 0.000882587, 0.002647761],
[ 0.005270463, 0.989459074, 0.002635231, 0.002635231, 0.000000000],
[ 0.000000000, 0.000000000, 0.005282977, 0.000000000, 0.000000000].
[ 0.000000000, 0.000000000, 0.000000000, 1.000000000, 0.000000000],
[ 0.000000000, 0.000000000, 0.000000000, 0.000000000, 1.000000000]
);

# It is simply done by assigning the value in @real_value given
# the corresponding (@row,@index) as coordinate;
# e.g. to assign \$real_value[-1] into
# \$actual_mat[\$row_index[-1]-1,\$col_index[-1]-1);
My question is: can we do matrix multiplication by just using sparse matrix representation above instead of @actual_mat.
Given the multiplicator
```my @p = (0.4,0.2,0.2,0.2,0.2);
we expect to get from p * SparseM
```my \$result = [
0.3989409,
0.2010541,
0.2000000,
0.2000000,
0.2000000,
];

---
neversaint and everlastingly indebted.......

Replies are listed 'Best First'.
Re: Sparse Matrix Multiplication Problem
by moritz (Cardinal) on Jan 27, 2009 at 11:42 UTC
If you're doing matrix multiplication yourself, you're doing it wrong. Use PDL::Sparse or another library for that.

(Usually sparse matrices are stored a bit different: for every row and every column there's a linked list of non-zero items together with their indexes, so that walking all non-zero elements is fast, ie following the links).

I'd like to quickly jump in to advocate Math::MatrixSparse, a pure-Perl sparse matrix library. It does just what its name promises. Note, however, that the original code has some annoying bugs, therefore, I maintain a separate version here.

Maybe you should ask the author for permissions to maintain the official module on CPAN? If you can't reach the author, have a look at the official PAUSE FAQ on this and mail the modules mailing list.

Cheers,
Steffen

Re: Sparse Matrix Multiplication Problem
by scorpio17 (Canon) on Jan 27, 2009 at 15:12 UTC

I usually reach for PDL in situations like this.

(Wow, it's like deja vu all over again...)

Re: Sparse Matrix Multiplication Problem
by gone2015 (Deacon) on Jan 27, 2009 at 17:17 UTC

So, you want to multiply your sparse array by a not-sparse vector @p. Seems straightforward enough:

```  my @r = (0) x scalar(@p) ;

for my \$i (0..\$#row_index) {
\$r[\$row_index[\$i]-1] += \$realval[\$i] * \$p[\$col_index[\$i]-1] ;
} ;

Dear Oshalla,
Just to let you know. I've ported your ingenious code to C++. They are so indispensable! Thanks a million again.
```#include <iostream>
#include <vector>
#include <boost/assign/std/vector.hpp>
#include <vector>
#include <stdio.h>
#include <stdlib.h>
using namespace boost::assign;
using namespace std;

int main  ( int arg_count, char *arg_vec[] ) {
vector <double> realValue;
realValue +=
0.994704478,
0.989459074,
0.994717023,
1.000000000,
1.000000000,
0.002647761,
0.005282977,
0.000882587,
0.005270463,
0.000882587,
0.002635231,
0.000882587,
0.002635231;

vector<int>rowIndex;
rowIndex +=  1, 2, 3, 4, 5, 1, 3, 1, 2, 1, 2, 1, 2;

vector<int>colIndex;
colIndex  +=   1, 2, 3, 4 ,5, 5, 3, 2, 1, 3, 3, 4, 4;

vector<int>matDim;
matDim  += 5,5; // M always equal to N

vector <double> theP;
theP += 0.4, 0.2, 0.2, 0.2, 0.2;

vector <double> Result;
Result.assign(theP.size(), 0);

for (int i= 0; i < rowIndex.size(); i++) {
Result[rowIndex[i]-1] += realValue[i] * theP[colIndex[i]-1];
}

// print it
for (int k =0; k < Result.size(); k++){
cout << Result[k] << endl;
}

return 0;
}

---
neversaint and everlastingly indebted.......
Dear oshalla,

I also tried to modify your code to do multiplication of transpose of the sparse matrix. With the following code:
```use Data::Dumper;
my @realval = (
0.994704478,
0.989459074,
0.994717023,
1.000000000,
1.000000000,
0.002647761,
0.005282977,
0.000882587,
0.005270463,
0.000882587,
0.002635231,
0.000882587,
0.002635231,);

my @row_index = qw( 1 2 3 4 5 1 3 1 2 1 2 1 2);
my @col_index =  qw(1 2 3 4 5 5 3 2 1 3 3 4 4);
my @mat_dim  = (5,5); # M always equal to N

my @trans_r = (0) x scalar(@p);

for my \$j (0..\$#row_index) {
\$trans_r[\$col_index[\$i]-1] += \$realval[\$i] * \$p[\$row_index[\$i]-1] ;
}

print Dumper \@trans_r;
Not sure why it doesn't give this desired result:
```@result = (0.3989359, 0.1982448, 0.2008801, 0.2008801, 0.2010591);
Truly need your expert advice on this. Thanks and hope to hear from you again.

---
neversaint and everlastingly indebted.......

You really should make friends with use strict and its cousin use warnings... With the code as given:

• use strict complains about \$i, which is used in the loop, but doesn't seem to be declared anywhere

-- that's what it means by 'Global symbol "\$i" requires explicit package name'... we just have to accept that when looked at from the correct angle, this makes perfect sense.

• use warnings would warn you that each time \$i is used it is uninitialised (though with use strict you wouldn't get that far !).</c>

Why not use a nested hash?
by Eckhart (Initiate) on Aug 27, 2009 at 13:18 UTC
Dear monks,

please accept my sincere apologies if you find the following question trivial.

Could anyone tell me why a self-made sparse matrix multiplication without PDL is bound to fail, as moritz said? In my opinion, anyone can use twice-nested hashs in order to represent a sparse 2d-matrix:

```\$matrix{\$row}{\$column}=\$value
This doesn't store the zeroes, and multiplication is also straightforward, by looping over appropiate key lists.

I don't see why this solution should be slower than using piddles.

Thank you very much for shedding some light on this issue.

Create A New User
Node Status?
node history
Node Type: perlquestion [id://739156]
Approved by citromatik
Front-paged by monkfan
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others having an uproarious good time at the Monastery: (9)
As of 2020-02-28 15:58 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
What numbers are you going to focus on primarily in 2020?

Results (124 votes). Check out past polls.

Notices?