Vectorize Matrix Multiplication

[Ben Desfosses]

I have many 3x3 matrices stored in the form: C(1:3,1:3,1:n)
where n is large.

I want to perform a cumulative multiplication of this set
of matrices.

If C was a vector C(1:n) I could just use: cumprod(C)

I tried cumprod(C, 3) and it gives me an element by element
cumulative product not matrix by matrix.

Can someone tell me how I might process or restructure this
set of data to get my desired result?



Thanks for any assistance.

[Peter Boettcher]

Post by Ben Desfosses
I have many 3x3 matrices stored in the form: C(1:3,1:3,1:n)
where n is large.
I want to perform a cumulative multiplication of this set
of matrices.
If C was a vector C(1:n) I could just use: cumprod(C)
I tried cumprod(C, 3) and it gives me an element by element
cumulative product not matrix by matrix.

My mex file ndfun does this.

This should do it:

D = ndfun('mprofd', C)


Find it at http://www.mit.edu/~pwb/matlab/ndfun. You'll have to compile
it yourself.


-Peter

[Ben Desfosses]

Post by Peter Boettcher

Post by Ben Desfosses
I have many 3x3 matrices stored in the form: C

(1:3,1:3,1:n)

Post by Peter Boettcher

Post by Ben Desfosses
where n is large.
I want to perform a cumulative multiplication of this

set

Post by Peter Boettcher

Post by Ben Desfosses
of matrices.
If C was a vector C(1:n) I could just use: cumprod(C)
I tried cumprod(C, 3) and it gives me an element by

element

Post by Peter Boettcher

Post by Ben Desfosses
cumulative product not matrix by matrix.

My mex file ndfun does this.
D = ndfun('mprofd', C)
Find it at http://www.mit.edu/~pwb/matlab/ndfun. You'll

have to compile

Post by Peter Boettcher
it yourself.
-Peter

Thanks so much, "ndfun" will be a life saver… if I can get
it to run

I compiled your c code into a .mexw32 and now I'm supposed
to be able to call it like any other .m file but I get an
error message saying, “Unknown command”, when I call it in
my .m file.

I had the ndfun.c file in my working directory.
I typed: mex ndfun.c C:\MATLAB\R2007a\extern\lib\win32
\lcc\libmwlapack.lib
It created: ndfun.mexw32 in my working directory
I tried to call: D = ndfun('mprofd', num2str(C));
And it returns: ??? Unknown command

Should I start a new thread for this problem???

Thanks,
-Ben

[Peter Boettcher]

Post by Ben Desfosses

Post by Peter Boettcher

I have many 3x3 matrices stored in the form: C (1:3,1:3,1:n) where
n is large.
I want to perform a cumulative multiplication of this set of
matrices.

My mex file ndfun does this.
D = ndfun('mprofd', C)

I had the ndfun.c file in my working directory.
I typed: mex ndfun.c C:\MATLAB\R2007a\extern\lib\win32
\lcc\libmwlapack.lib
It created: ndfun.mexw32 in my working directory
I tried to call: D = ndfun('mprofd', num2str(C));
And it returns: ??? Unknown command

Sorry, I made a typo. "Unknown command" is a message from NDFUN,
complaining about the string 'mprodf'. 'mprod' is the correct command.
You should also download ndfun.m from the same location (put it in the
same directory). It contains the help text, so that "help ndfun" will
give you the full usage information.

-Peter

[Ben Desfosses]

I should have caught that myself,

Thanks

[Ben Desfosses]

Nothing ever goes smoothly…

Peter,

Before I try I just wanted to get your opinion.
Do you think it will be simple enough if I try to modify
your C program to calculate A(:,:,2)*A(:,:,1) instead of A
(:,:,1)*A(:,:,2)? And also return all the calculated
matrices and not just the final result?

I’m not a big C programmer so I didn’t want to dive into a
futile task before asking.

The for loop for this is fast enough if I don't store each
Matrix but when I try to keep them the time increases by a
huge amount.

Thanks,
Ben

[Peter Boettcher]

Post by Ben Desfosses
Nothing ever goes smoothly…
Peter,
Before I try I just wanted to get your opinion.
Do you think it will be simple enough if I try to modify
your C program to calculate A(:,:,2)*A(:,:,1) instead of A
(:,:,1)*A(:,:,2)? And also return all the calculated
matrices and not just the final result?

Just flip your matrix before running the program. But as you suggest,
changing it in C code wouldn't be bad. Line 402 is the start of the
small loop that runs over the third dimension. It starts at 1 and goes
to n. Reversing that loop should do it.

Post by Ben Desfosses
I’m not a big C programmer so I didn’t want to dive into a
futile task before asking.

Storing all the matrices would be a much bigger modification. You could
certainly do it, but if you don't know much C it might be out of reach.

The current code does the following, where numbers represent the matrix
at that page:

1*2 -> dst
dst*3 -> tmp
tmp*4 -> dst
dst*5 -> tmp
...

The multiply routine can't have the same input and output space, which
is why we need to flip back and forth between two output areas. To do
you what you're looking for, you'd need to allocate a bigger output
array, and fix the indexing. Then, where s is source and d is destination,

1s->1d (just copy)
1d*2s-> 2d
2d*3s-> 3d
...

Post by Ben Desfosses
The for loop for this is fast enough if I don't store each
Matrix but when I try to keep them the time increases by a
huge amount.

Show us the "for loop". There might be an enhancement to help at that
level.

-Peter

[Ben]

n=size(A, 3);


Yor Program does:
C = A(:,:,1);
for i=2:n
C = C*A(:,:,i);
end


The first change would be to change to this:
C = A(:,:,1);
for i=2:n
C = A(:,:,i)*C;
end

This is What I realy need…
C(:,:,1)=A(:,:,1);
for i=2:n
C(:,:,i)=A(:,:,i)*C(:,:,i-1);
end

Based on you comment about the loop I looked at it my self
and realized that the “pre-allocate for speed” auto comment
actually makes a difference and changed it to:

C=A;
for i=2:n
C(:,:,i)=A(:,:,i)*C(:,:,i-1);
end

This made a huge difference. Rookie mistake
Thanks for the help.

[Ben Desfosses]

Correction

I called:
D = ndfun('mprofd', C);

not
D = ndfun('mprofd', num2str(C));

I was trying somthing I saw in a different post to correct
the problem

[Kenny]

Here's code that doesn't require outside libraries:


% MATMULT Vectorized matrix multiplication.
%
% Usage: C = matmult(A, B)
% A is an l x m x k list of k matrices
% B is an m x n x k list of k matrices or an m x k list of k vectors
%
function C = matmult(A, B)
if length(size(B))==2
C = matmult_helper(A, B);
else
C = cell2mat(arrayfun( @(k) reshape(matmult_helper(A, squeeze(B(:,k,:))),[size(A,1),1,size(B,3)]), 1:size(B,2), 'UniformOutput', false));
end;


function y = matmult_helper(A, x)

Post by Ben Desfosses
I have many 3x3 matrices stored in the form: C(1:3,1:3,1:n)
where n is large.
I want to perform a cumulative multiplication of this set
of matrices.
If C was a vector C(1:n) I could just use: cumprod(C)
I tried cumprod(C, 3) and it gives me an element by element
cumulative product not matrix by matrix.
Can someone tell me how I might process or restructure this
set of data to get my desired result?
Thanks for any assistance.

[Matt J]

That's an interesting implementation if you are looking to duplicate some of the functionality of MTIMESX

http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

and is even quite competitive with it for small matrices. It's not what the OP asked for however. The OP is looking for a cumulative matrix product over the slices
A(:,:,k), so there shouldn't even be two arguments A and B.

It might also be worth clarifying for posterity readers that matmult is not expected to outperform even a basic for-loop when size(A,2) gets sufficiently large. The following comparisons confirm this


A=rand(50,50,200);

tic; C0=forloopBaseline(A,A); toc %Elapsed time is 0.039035 seconds.

tic; C1=matmult(A,A); toc %Elapsed time is 0.292982 seconds.

tic; C2=mtimesx(A,A); toc %Elapsed time is 0.008018 seconds.



function C=forloopBaseline(A,B)

n=size(A,3);

C=A;
for i=1:n
C(:,:,i)=A(:,:,i)*B(:,:,i);
end



function C = matmult(A, B)
if length(size(B))==2
C = matmult_helper(A, B);
else
C = cell2mat(arrayfun( @(k) reshape(matmult_helper(A, squeeze(B(:,k,:))),[size(A,1),1,size(B,3)]), 1:size(B,2), 'UniformOutput', false));
end;


function y = matmult_helper(A, x)
y = squeeze(sum(bsxfun(@times, A, reshape(x, [1,size(A,2),size(x,2)])),2));

[James Tursa]

Post by Matt J
That's an interesting implementation if you are looking to duplicate some of the functionality of MTIMESX
http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support
and is even quite competitive with it for small matrices. It's not what the OP asked for however. The OP is looking for a cumulative matrix product over the slices
A(:,:,k), so there shouldn't even be two arguments A and B.

FYI, the next release of MTIMESX will include the matrix multiply equivalents of prod and cumprod for MxMxK arrays. The prototype is working, but I haven't found the time yet to create appropriate test suites so haven't uploaded it to the FEX.

James Tursa

[Matt J]

Post by Kenny
% MATMULT Vectorized matrix multiplication.
%
% Usage: C = matmult(A, B)
% A is an l x m x k list of k matrices
% B is an m x n x k list of k matrices or an m x k list of k vectors

======================

Also, even for small size(B,2), you can do a bit better by implementing the same technique used by MATMULT with a for-loop over B(:,k,:)

A=rand(5,5,1e6);

tic; C1=matmult(A,A); toc %Elapsed time is 2.894654 seconds.
tic; C2=kforloop(A,A); toc %Elapsed time is 2.763034 seconds.


function C=kforloop(A,B)

n=size(A,2);
m=size(A,3);
C=A;

for k=1:n
C(:,k,:)=sum(bsxfun(@times,A,reshape(B(:,k,:),1,n,m)),2);
end