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

Sorry if this has already been asked before, but I searched and could not find anything... I was trying to multiply two fairly large matrices (~5000 by 5000) using the MatrixReal module. It ended up being so slow that I actually decided to write the matrices to a file, use a Fortran code to read, multiply and write the result, and finally read the result back. Is there an elegant solution to my problem??? Thanks a lot, Ludovic

Replies are listed 'Best First'.
Re: Fast matrix multiplication
by hardburn (Abbot) on Aug 11, 2003 at 19:10 UTC

    Interestingly, this came up today on perl.com's entry from the Perl Cookbook.

    Use PDL, which overloads the '*' operator (among others) to multiply a matrix.

    ----
    I wanted to explore how Perl's closures can be manipulated, and ended up creating an object system by accident.
    -- Schemer

    Note: All code is untested, unless otherwise stated

[reply]
      In PDL, '*' is for scalar multiplication, which automatically "broadcasts" up so that two ndarrays can be scalar-multiplied so the result is each of the six "slots" being the product of the two numbers. Matrix multiplication is the overload of the 'x' operator. 
      pdl> p $x = sequence(2,2)
[
 [0 1]
 [2 3]
]

pdl> p $y = sequence(2,2) + 3
[
 [3 4]
 [5 6]
]

pdl> p $x * $y
[
 [ 0  4]
 [10 18]
]

pdl> p $x x $y
[
 [ 5  6]
 [21 26]
]

      [download]
[reply]
[d/l]
Re: Fast matrix multiplication
by Zaxo (Archbishop) on Aug 11, 2003 at 19:39 UTC

    A 5000 by 5000 matrix contains 25 million elements. That is 200 million bytes of 64 bit doubles, together with 5001 perl array header structures. You need three of those monsters to do an array multiplication. Do you suppose you're hitting swap?

    Since most N-element matrix multiplication schemes involve O(N**3) multiplications, you're on the bad side of a scaling problem, as well.

    Your strategy of serializing to disk looks pretty good in light of that. You may want to transpose the right hand matrix in order to get the columns as a single array ref.

    The best cure would be less brute force. Try to reduce your problem to a block diagonal form if you can. Think about transformations that may simplify things. Would approximate knowledge of the eigenvectors help?

    In any case, PDL is superior to pure perl math packages for speed and flexibility.

    After Compline,
    Zaxo

[reply]
Re: Fast matrix multiplication
by dragonchild (Archbishop) on Aug 11, 2003 at 19:02 UTC
    That is an elegant solution. There is nothing that says you cannot use Fortran to do something Fortran is good at and Perl isn't.

    Now, I would launch the Fortran executable from within the Perl, but that's just me. :-)

    ------
    We are the carpenters and bricklayers of the Information Age.

    The idea is a little like C++ templates, except not quite so brain-meltingly complicated. -- TheDamian, Exegesis 6

    Please remember that I'm crufty and crochety. All opinions are purely mine and all code is untested, unless otherwise specified.

[reply]
Re: Fast matrix multiplication
by blokhead (Monsignor) on Aug 11, 2003 at 19:13 UTC
    This may be a good time to learn Inline::C (is there an Inline::Fortran?). It will save you the time of writing a large matrix to a file and reloading it in another program, etc. You may also want to take a look at PDL, as it claims to be very fast for manipulating large matrix. I've never used it, but the interface looks pretty logical, and it may be simpler than to try to roll your own super-fast matrix mangler in C.

    blokhead

[reply]
Re: Fast matrix multiplication
by Abigail-II (Bishop) on Aug 11, 2003 at 21:42 UTC
    Perl is pretty bad when it comes to doing large scale arithmatic, and Perl is really lousy when it comes to storing millions of scalars. If using Fortran works for you, don't change it. After all, Perl is a glue language and I've never understood those people who always strive to have a pure Perl solution. But if you want to look further than your use of Fortran, you should look into PDL.

    Abigail

[reply]
Re: Fast matrix multiplication
by bean (Monk) on Aug 11, 2003 at 21:02 UTC
    If you have a sparse matrix (mostly zeros) there are algorithms out there that will speed things up considerably. Google "sparse matrix multiplication". If not, you'll probably have to hook in some compiled code (as you have done) - to get a little extra speed you might try MMX extensions.
    Update:
    Here's the basic idea behind sparse matrix multiplication. Collapse the vectors into lists of the non-zero elements and corresponding lists of their indexes in the original vectors. Then traverse the lists of indexes for the vectors you are multiplying - you only have to multiply the vector elements when their indexes match up. Your result is the intersection of the lists of indexes and the multiplied elements of those indexes. (Obviously, this only helps if most of the elements are zeros - however, that's fairly common in the sciences.)
[reply]
Re: Fast matrix multiplication
by waswas-fng (Curate) on Aug 11, 2003 at 20:02 UTC
    you could try benchmarking Math::Cephes::Matrix. It is a perl frontend to the Cephes optimized math lib (written in C) I have found it to be very fast.

    -Waswas
[reply]
Re: Fast matrix multiplication
by mattr (Curate) on Aug 13, 2003 at 11:34 UTC
    I have not used this, and I do not know about the raw bandwidth issue, but you may like to try GSL (Gnu Scientific Library) with ATLAS for enhanced BLAS matrix operations. See "Alternative BLAS Libraries" on the GSL page. Perhaps people who do more matrix multiplication can compare this to the other two packages mentioned above.

    However the GSL page mentions PDL::GSL which I can't seem to find. And I don't know if Math::GSL covers matrix multiplication (you'd think it would..). On the other hand, if the important part is the BLAS implementation and you are using ATLAS, why not just use ATLAS either with Inline::C or better yet from command line, skipping all the to-from perl overhead.

    From http://math-atlas.sourceforge.net/,
    The ATLAS (Automatically Tuned Linear Algebra Software) project is an ongoing research effort focusing on applying empirical techniques in order to provide portable performance. At present, it provides C and Fortran77 interfaces to a portably efficient BLAS implementation, as well as a few routines from LAPACK.

    If someone else could comment on this I would appreciate it. Very likely the PDL people are well aware of ATLAS too.

[reply]
      PDL has bindings to parts of GSL, rather than a single big PDL::GSL. See https://metacpan.org/dist/PDL and text-search for "gsl". There are indeed plans to use BLAS, or alternative implementations such as ATLAS, more pervasively.
[reply]
      P.S. The ATLAS homepage says it is used as an option in R.
[reply]
Re: Fast matrix multiplication
by mattr (Curate) on Aug 16, 2003 at 13:45 UTC
    Incidentally I just saw some information about using graphics processors on 3d accelerated video cards for simulation including matrix calculations.. seems to be coming of age. You may wish to check out:

    paper at AT&T
    Torsten Moller

    the post (on slashdot) was by Amit Patel.

[reply]
Re: Fast matrix multiplication
by Anonymous Monk on Aug 12, 2003 at 15:02 UTC
    If you plan on doing a lot of matrix work consider shelling out to R. R is especially designed for this, and handles complex matrix operations fairly efficiently and without too much coding.
[reply]

Back to Seekers of Perl Wisdom