Beefy Boxes and Bandwidth Generously Provided by pair Networks
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 ( #739156=perlquestion: print w/replies, xml ) 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,
      Thank you so much for your reply. It is truly useful.

      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.

Log In?
Username:
Password:

What's my password?
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?