Thread Subject: Is there a parallel version of matrix multiplication in MATLAB?

Hi, I am calculating powers of large matrix. For example, taking a matrix of size 2000*2000 to the power of 6. I know that MATLAB is known for its performance in matrix operations. But in this case, it is still quite slow.

I have the parallel toolbox, I'm wondering is there any parallel version of matrix multiplication available? If so, could anyone give me a short example of how to use it?

Thanks so much!
Cindy

"Cindy" <xzhan2@fhcrc.org> wrote in message <hcqehk$17u$1@fred.mathworks.com>...
> Hi, I am calculating powers of large matrix. For example, taking a matrix of size 2000*2000 to the power of 6. I know that MATLAB is known for its performance in matrix operations. But in this case, it is still quite slow.
>
> I have the parallel toolbox, I'm wondering is there any parallel version of matrix multiplication available? If so, could anyone give me a short example of how to use it?
>
> Thanks so much!
> Cindy

MATLAB uses calls to a BLAS library to do matrix multiplication, and those calls use parallel processing. I believe this became part of MATLAB beginning in R2007a (although you had to manually turn it on in that version), and became the default behavior in R2008a(?). So you will not get any speedup by trying to manually piece up the calculation and parallelize (is that a word?) it yourself.

James Tursa

"Cindy" <xzhan2@fhcrc.org> writes:

> Hi, I am calculating powers of large matrix. For example, taking a matrix of
> size 2000*2000 to the power of 6. I know that MATLAB is known for its
> performance in matrix operations. But in this case, it is still quite slow.
>
> I have the parallel toolbox, I'm wondering is there any parallel version of
> matrix multiplication available? If so, could anyone give me a short example
> of how to use it?

There's matrix multiplication defined for distributed arrays that ship with the
Parallel Computing Toolbox. On a single machine (i.e. using "matlabpool local"),
the multithreaded version of mtimes is probably faster, but here's how to do it
(using R2009b syntax):

matlabpool open % plus any options
spmd
  d = codistributed.rand( 2000 );
  d * d * d * d * d * d;
end

Cheers,

Edric.

"Cindy" <xzhan2@fhcrc.org> wrote in message
news:hcqehk$17u$1@fred.mathworks.com...
> Hi, I am calculating powers of large matrix. For example, taking a matrix
> of size 2000*2000 to the power of 6. I know that MATLAB is known for its
> performance in matrix operations. But in this case, it is still quite
> slow.

From the way you phrased your question, I suspect you're raising the matrix
to the 6th power using repeated multiplication. If so, why? Just use ^
(the MPOWER operator) instead.

a = rand(2000);
tic
    b = a*a*a*a*a*a;
t1 = toc
tic
    c = a^6;
t2 = toc

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

"Steven Lord" <slord@mathworks.com> wrote in message <hcs5it$2gu$1@fred.mathworks.com>...
>
> "Cindy" <xzhan2@fhcrc.org> wrote in message
> news:hcqehk$17u$1@fred.mathworks.com...
> > Hi, I am calculating powers of large matrix. For example, taking a matrix
> > of size 2000*2000 to the power of 6. I know that MATLAB is known for its
> > performance in matrix operations. But in this case, it is still quite
> > slow.
>
> From the way you phrased your question, I suspect you're raising the matrix
> to the 6th power using repeated multiplication. If so, why? Just use ^
> (the MPOWER operator) instead.
>
> a = rand(2000);
> tic
> b = a*a*a*a*a*a;
> t1 = toc
> tic
> c = a^6;
> t2 = toc
>
> --

One could use the binary decomposition of the power to perform computation with less calculation:

a = rand(2000);
n = 6;

% Compute a^n

% straightforward multiplication
tic
    b = a;
    for k=2:n
        b = a*b;
    end
t1 = toc % 3.6662 seconds

% Matlab power
tic
    c = a^n;
t2 = toc % 3.0213 seconds

% binary decomposition
tic
n = 6;
nb = dec2bin(n);
ak = a;
d = 1;
for k=length(nb)-1:-1:1
    ak = ak*ak;
    if nb(k);
        d = d*ak;
    end
end
t3 = toc % 2.2377 seconds

% Bruno

% Sorry fr the bug

tic
n = 6;
nb = dec2bin(n)-'0';
ak = a;
d = 1;
for k=length(nb)-1:-1:1
    ak = ak*ak;
    if nb(k)
        d = d*ak;
    end
end
t3 = toc % 2.2377 seconds

% Bruno

Arrg, there is still another bug...

% binary decomposition
nb = dec2bin(n)-'0';
ak = a;
if nb(end)
    d = a;
else
    d = 1;
end
for k=length(nb)-1:-1:1
    ak = ak*ak;
    if nb(k)
        d = d*ak;
    end
end


Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hcs9ef$c0u$1@fred.mathworks.com>...
> Arrg, there is still another bug...
>
> % binary decomposition
> nb = dec2bin(n)-'0';
> ak = a;
> if nb(end)
> d = a;
> else
> d = 1;
> end
> for k=length(nb)-1:-1:1
> ak = ak*ak;
> if nb(k)
> d = d*ak;
> end
> end
>
>
> Bruno

Hey Bruno ... out of curiosity, do you happen to know if the polyvalm function uses this type of scheme (or similar) to minimize total computations? I could code something up and make some timing tests, but wondered if you already knew the answer.

James Tursa

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hct1kq$kd$1@fred.mathworks.com>...

>
> Hey Bruno ... out of curiosity, do you happen to know if the polyvalm function uses this type of scheme (or similar) to minimize total computations? I could code something up and make some timing tests, but wondered if you already knew the answer.

Hi James,

I believe (!) POLYVAL uses horner's scheme. I have no idea how they compare to each other. I would say Horner's is more efficient in general, but if you confirm that it will be great!!! Thanks for bringing this interesting question up.

Bruno

Sorry for the typo:

POLYVAL -> POLYVALM

Bruno

"Cindy" <xzhan2@fhcrc.org> wrote in message <hcqehk$17u$1@fred.mathworks.com>...
> Hi, I am calculating powers of large matrix. For example, taking a matrix of size 2000*2000 to the power of 6. I know that MATLAB is known for its performance in matrix operations. But in this case, it is still quite slow.
>
> I have the parallel toolbox, I'm wondering is there any parallel version of matrix multiplication available? If so, could anyone give me a short example of how to use it?
>
> Thanks so much!
> Cindy



Cathy,

The answer to your question depends on what you want to do with A^6. If you must calculate A^6 then Stephen Lord points out that A^6 is faster than A*A*A*A*A*A, because A^6 uses repeated squaring: A2 = A*A, A4 = A2*A2, A6 = A2*A4, which use 3 mat-mults instead of 5.

However, if you want to solve A^k = b then Golub & Van Loan's method (page 121) is much faster:

[L,U,P] = lu(A);
for m = 1:k
    y = L\(P*b);
    b = y;
    x3 = U\y;
    b = x3;
end;

This has complexity O(n^3 + k*n^2), which is much faster than Matlab's O(n^3*log_2 k), and yours O(n^3*k).

See http://www.derekroconnor.net/NA/LE/LE-2006-6.pdf

Derek O'Connor

"Derek O'Connor" <derekroconnor@eircom.net> wrote in message <hcun0e$ngj$1@fred.mathworks.com>...
> "Cindy" <xzhan2@fhcrc.org> wrote in message <hcqehk$17u$1@fred.mathworks.com>...
> > Hi, I am calculating powers of large matrix. For example, taking a matrix of size 2000*2000 to the power of 6. I know that MATLAB is known for its performance in matrix operations. But in this case, it is still quite slow.
> >
> > I have the parallel toolbox, I'm wondering is there any parallel version of matrix multiplication available? If so, could anyone give me a short example of how to use it?
> >
> > Thanks so much!
> > Cindy
>
>
>
> Cathy,
>
> The answer to your question depends on what you want to do with A^6. If you must calculate A^6 then Stephen Lord points out that A^6 is faster than A*A*A*A*A*A, because A^6 uses repeated squaring: A2 = A*A, A4 = A2*A2, A6
> snip

Sorry that should be Cindy, not Cathy.

 Derek O'Connor

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hct2k7$ao$1@fred.mathworks.com>...
> "James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hct1kq$kd$1@fred.mathworks.com>...
>
> >
> > Hey Bruno ... out of curiosity, do you happen to know if the polyvalm function uses this type of scheme (or similar) to minimize total computations? I could code something up and make some timing tests, but wondered if you already knew the answer.
>
> Hi James,
>
> I believe (!) POLYVAL uses horner's scheme. I have no idea how they compare to each other. I would say Horner's is more efficient in general, but if you confirm that it will be great!!! Thanks for bringing this interesting question up.
>
> Bruno

Well, Horner's scheme is just a simple nested rewrite of the polynomial evaluation. That is definitely not the most computationally efficient method if a lot of the polynomial coefficients are zero. e.g., for A^16 Horner's scheme, if taken literally, would do 15 multiplies, whereas the binary decomposition scheme would do only 4 multiplies. But for a general case where there is a mix of zero and non-zero coefficients, it seems you could in some cases still use some form of binary decomposition to get the minimum overall computation, but at the expense of storing more intermediate matrices. I may have to make some tests ...

James Tursa

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message
news:hct1kq$kd$1@fred.mathworks.com...

*snip*

> Hey Bruno ... out of curiosity, do you happen to know if the polyvalm
> function uses this type of scheme (or similar) to minimize total
> computations? I could code something up and make some timing tests, but
> wondered if you already knew the answer.

Since POLYVAL and POLYVALM (depending on which you're interested in) are
both M-files, you can simply read the code and check. POLYVALM is fairly
straightforward; POLYVAL is a little more complicated because it needs to
handle cases that POLYVALM doesn't (like the MU and S outputs from POLYFIT
and the DELTA output it returns), but it's still not that bad.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

"Steven Lord" <slord@mathworks.com> wrote in message <hcuqqe$q72$1@fred.mathworks.com>...
>
> Since POLYVAL and POLYVALM (depending on which you're interested in) are
> both M-files, you can simply read the code and check.

Ah ... Horner's method.

James Tursa

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hcuqfj$4m7$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hct2k7$ao$1@fred.mathworks.com>...
> > "James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hct1kq$kd$1@fred.mathworks.com>...
> >
> > >
> > > Hey Bruno ... out of curiosity, do you happen to know if the polyvalm function uses this type of scheme (or similar) to minimize total computations? I could code something up and make some timing tests, but wondered if you already knew the answer.
> >
> > Hi James,
> >
> > I believe (!) POLYVAL uses horner's scheme. I have no idea how they compare to each other. I would say Horner's is more efficient in general, but if you confirm that it will be great!!! Thanks for bringing this interesting question up.
> >
> > Bruno
>
> Well, Horner's scheme is just a simple nested rewrite of the polynomial evaluation. That is definitely not the most computationally efficient method if a lot of the polynomial coefficients are zero. e.g., for A^16 Horner's scheme, if taken literally, would do 15 multiplies, whereas the binary decomposition scheme would do only 4 multiplies. But for a general case where there is a mix of zero and non-zero coefficients, it seems you could in some cases still use some form of binary decomposition to get the minimum overall computation, but at the expense of storing more intermediate matrices. I may have to make some tests ...

Yup, but in case of a polynomial has all non-zeros coefficients, I believe the number of matrix multiplications needed is the same for both Horner's and binary-accumulation.

But if you make some test, I would be interested in knowing the result.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hcuvi6$2m5$1@fred.mathworks.com>...
> "James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hcuqfj$4m7$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hct2k7$ao$1@fred.mathworks.com>...
> > > "James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hct1kq$kd$1@fred.mathworks.com>...
> > >
> > > >
> > > > Hey Bruno ... out of curiosity, do you happen to know if the polyvalm function uses this type of scheme (or similar) to minimize total computations? I could code something up and make some timing tests, but wondered if you already knew the answer.
> > >
> > > Hi James,
> > >
> > > I believe (!) POLYVAL uses horner's scheme. I have no idea how they compare to each other. I would say Horner's is more efficient in general, but if you confirm that it will be great!!! Thanks for bringing this interesting question up.
> > >
> > > Bruno
> >
> > Well, Horner's scheme is just a simple nested rewrite of the polynomial evaluation. That is definitely not the most computationally efficient method if a lot of the polynomial coefficients are zero. e.g., for A^16 Horner's scheme, if taken literally, would do 15 multiplies, whereas the binary decomposition scheme would do only 4 multiplies. But for a general case where there is a mix of zero and non-zero coefficients, it seems you could in some cases still use some form of binary decomposition to get the minimum overall computation, but at the expense of storing more intermediate matrices. I may have to make some tests ...
>
> Yup, but in case of a polynomial has all non-zeros coefficients, I believe the number of matrix multiplications needed is the same for both Horner's and binary-accumulation.
>
> But if you make some test, I would be interested in knowing the result.
>
> Bruno

I was curious, so I coded up a version of polyvalm that used a binary decomosition scheme for the powers. As expected, when the coefficients are all non-zero it ran a bit slower (more overhead but still had to do the same number of multiplies). At the other end, where there were a lot of zero coefficients, the savings showed up. In the extreme case where only the first coefficient was non-zero, my polyvalm2 function ran in less than 1/2 the time of polyval for a large (1500x1500) matrix, essentially matching the performance of mpower in this case. I don't know where the line is (number of zero coefficients) for determining which method is faster. Certainly Horner's scheme is more memory efficient, so it is the more robust method. Maybe I will clean it up my code and post it to the FEX.

James Tursa

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hcv44h$lct$1@fred.mathworks.com>...

> I was curious, so I coded up a version of polyvalm that used a binary decomosition scheme for the powers. As expected, when the coefficients are all non-zero it ran a bit slower (more overhead but still had to do the same number of multiplies). At the other end, where there were a lot of zero coefficients, the savings showed up. In the extreme case where only the first coefficient was non-zero, my polyvalm2 function ran in less than 1/2 the time of polyval for a large (1500x1500) matrix, essentially matching the performance of mpower in this case. I don't know where the line is (number of zero coefficients) for determining which method is faster. Certainly Horner's scheme is more memory efficient, so it is the more robust method. Maybe I will clean it up my code and post it to the FEX.
>
> James Tursa

Thanks James. The conclusion of your test meet indeed our expectation. Good point about the stability. Personally, I would not use the binary-power scheme accept for the mono term polynomial for this very reason. Still it's good to have an alternative algorithm available.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hcv52f$kbv$1@fred.mathworks.com>...
> "James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hcv44h$lct$1@fred.mathworks.com>...
>
> > I was curious, so I coded up a version of polyvalm that used a binary decomosition scheme for the powers. As expected, when the coefficients are all non-zero it ran a bit slower (more overhead but still had to do the same number of multiplies). At the other end, where there were a lot of zero coefficients, the savings showed up. In the extreme case where only the first coefficient was non-zero, my polyvalm2 function ran in less than 1/2 the time of polyval for a large (1500x1500) matrix, essentially matching the performance of mpower in this case. I don't know where the line is (number of zero coefficients) for determining which method is faster. Certainly Horner's scheme is more memory efficient, so it is the more robust method. Maybe I will clean it up my code and post it to the FEX.
> >
> > James Tursa
>
> Thanks James. The conclusion of your test meet indeed our expectation. Good point about the stability. Personally, I would not use the binary-power scheme accept for the mono term polynomial for this very reason. Still it's good to have an alternative algorithm available.
>
> Bruno

I need to modify my earlier statement. After examining the results a bit more closely, my polyval2 function is running about 20% faster than mpower for the extreme case. This was unexpected to me since I thought I was using essentially the same algorithm as mpower for this case. Puzzling. I will try to clean up the code & post to the FEX tonight.

James Tursa

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hcv8bf$glg$1@fred.mathworks.com>...

> I need to modify my earlier statement. After examining the results a bit more closely, my polyval2 function is running about 20% faster than mpower for the extreme case. This was unexpected to me since I thought I was using essentially the same algorithm as mpower for this case. Puzzling.

Well, the timing I posted yesterday support that. When I post I have not read the help of mpower, until Derek brought this up. I's unclear to me why mpower is slighly slower (about 30% on my timing).

Bruno

"Derek O'Connor" <derekroconnor@eircom.net> wrote in message <hcun0e$ngj$1@fred.mathworks.com>...
 ----- snip ---
> However, if you want to solve A^k = b then Golub & Van Loan's method (page 121) is much faster:
>
> [L,U,P] = lu(A);
> for m = 1:k
> y = L\(P*b);
> b = y;
> x3 = U\y;
> b = x3;
> end;
 

The above should be

 However, if you want to solve (A^k)x = b then Golub & Van Loan's method (page 121) is much faster:
 
 [L,U,P] = lu(A);
 for m = 1:k
     b = L\(P*b);
     x = U\b;
     b = x;
 end;

Derek O'Connor

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hcvbtp$1ns$1@fred.mathworks.com>...
> "James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hcv8bf$glg$1@fred.mathworks.com>...
>
> > I need to modify my earlier statement. After examining the results a bit more closely, my polyval2 function is running about 20% faster than mpower for the extreme case. This was unexpected to me since I thought I was using essentially the same algorithm as mpower for this case. Puzzling.
>
> Well, the timing I posted yesterday support that. When I post I have not read the help of mpower, until Derek brought this up. I's unclear to me why mpower is slighly slower (about 30% on my timing).
>
> Bruno

Bruno,

FYI, I have posted the files here in case you are interested:

http://www.mathworks.com/matlabcentral/fileexchange/25780-polyvalm2-a-faster-matrix-polynomial-evaluator

http://www.mathworks.com/matlabcentral/fileexchange/25782-mpower2-a-faster-matrix-power-function

For polyvalm I was able to examine the m-file and determine the basic sources of the speedup, which are documented in the help. After putting in the "look-ahead" feature to save processing, I have found that polyvalm2 always beats polyvalm, even if the P vector has all non-zeros. Also, for the case where the X matrix is real and P is complex, polyvalm basically builds up the result using complex matrix multiplies at each step, whereas polyvalm2 builds up the matrix powers using only real matrix multiplies and then applies the complex scalar multiply at the end. This alone easily cuts the time in 1/2.

For mpower there was no m-file to examine that contained the algorithm, so I had to guess. But a very telling test is as follows:

>> A = rand(2000);
>> tic;A^1;toc
Elapsed time is 3.075316 seconds.
>> tic;mpower2(A,1);toc
Elapsed time is 0.061539 seconds.

It appears that MATLAB's mpower *probably* starts with the identity matrix and then applies the various matrix powers. i.e., there is an unnecessary matrix multiply at the start of the algorithm. My guess is that is what is causing the slower timing of MATLAB's mpower vs my posted mpower2 function (using the binary decomposition of the power as discussed earlier in this thread).

James Tursa

Great! Thanks James. I'll check out your function later on and play with it.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hd9r8s$o91$1@fred.mathworks.com>...
> Great! Thanks James. I'll check out your function later on and play with it.
>
> Bruno

Bruno,

I just made another submission that you might be interested in, a mex routine that does matrix multiplication using BLAS calls and other custom code. I know, I know, MATLAB already does this. In fact it is highly optimized to do matrix multiplication as fast as possible, right? But I discovered that MATLAB does not always take advantage of some of the special cases (symmetric, conjugate, transpose, etc). My mtimesx code can beat MATLAB for speed in quite a few of these cases ... and not just by a few percentage points. I am talking about 2x - 3x (or more) improvement in speed. Results are highly dependent on computer, compiler, and MATLAB version. You can find it here:

http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-using-blas-calls

James Tursa

Great James,

I'll play with it and let you know how it performs under the various platforms available to me.

Later and cheers,

Bruno