Vectorization

[Junior]

Dear all,
I would like to vectorize the following piece of code

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

Any suggestion on how to go about that?

Thanks,
Junior

[Peter Boettcher]

"Junior " <junio...@gmail.com> writes:

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

Stick with the "for" loop, or try my MEX file ndfun from

http://www.mit.edu/~pwb/matlab/

A = ndfun('mult', B, C);

-Peter

[Matt Fig]

That for loop can be sped up by allowing the index of each loop to decrement instead of increase, assuming you didn't pre-allocate A.

for i=n:-1:1
for j=m:-1:1

[Junior]

"Matt Fig" <spam...@yahoo.com> wrote in message <gq8uc5$7n7$1...@fred.mathworks.com>...

Thanks for the answers guys.
MATT: do you mean I should not preallocate A?
The Mex solution is probably interesting as well but I would like to do it all in matlab directly.

Junior

[Bruno Luong]

Here is a fast version of 100% Matlab (no mex and no for-loop):

m = 3;
n = 4;
p = 2;
nb = 1e5;

% For loop
tic
A = rand(m,n,nb);
B = rand(n,p,nb);
C=zeros(m,p,nb);
for k=1:nb
C(:,:,k) = A(:,:,k) * B(:,:,k);
end
toc % 1.345668 seconds.

% Sparse matrix
tic
I = repmat(reshape(1:m*nb,m,1,nb),[1 n 1]);
J = repmat(reshape(1:n*nb,1,n,nb),[m 1 1]);
As = sparse(I(:),J(:),A(:));
I = repmat(reshape(1:n*nb,n,1,nb),[1 p 1]);
J = repmat(reshape(1:p*nb,1,p,nb),[n 1 1]);
Bs = sparse(I(:), J(:), B(:));
Cs = reshape(nonzeros(As * Bs),[m p nb]);
toc % 1.345668 seconds.

isequal(C, Cs) % <- Check

[Bruno Luong]

Sorry, respective timings are:
- For loop: Elapsed time is 1.342406 seconds.
- Sparse: Elapsed time is 0.332618 seconds.

Bruno

[Matt Fig]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <gq95ag$c82$1...@fred.mathworks.com>...

Since you use A and B in both timings, it seems unfair to compare their creation in the for loop timing. Even with them taken out, I still get that the sparse solution is faster.

% Create all data.

m = 3;
n = 4;
p = 2;
nb = 1e5;

A = rand(m,n,nb);
B = rand(n,p,nb);

% For loop
tic

C=zeros(m,p,nb);
for k=1:nb
C(:,:,k) = A(:,:,k) * B(:,:,k);
end

toc %Elapsed time is 0.409208 seconds.

% Sparse matrix
tic
I = repmat(reshape(1:m*nb,m,1,nb),[1 n 1]);
J = repmat(reshape(1:n*nb,1,n,nb),[m 1 1]);
As = sparse(I(:),J(:),A(:));
I = repmat(reshape(1:n*nb,n,1,nb),[1 p 1]);
J = repmat(reshape(1:p*nb,1,p,nb),[n 1 1]);
Bs = sparse(I(:), J(:), B(:));
Cs = reshape(nonzeros(As * Bs),[m p nb]);

toc %Elapsed time is 0.360580 seconds.

To the OP: No, I wasn't suggesting you not pre-allocate.

[Jinhui Bai]

"Junior " <junio...@gmail.com> wrote in message <gq8ism$rv$1...@fred.mathworks.com>...

Junior,
You might want to take a look at the discussions in another related post, where I refer to another solution based on "bsxfun(@times, ...)" and "sum." In other words, it is based on the direct definition of matrix product as pairwise inner products, although without additional memory usage due to the use of bsxfun.
http://www.mathworks.com/matlabcentral/newsreader/view_thread/244511#628393

In general, if the first two dimensions of A and B are small and/or the page dimension is large, these alternatives may speed things up. But if A and B are large and/or page dimension is small, it is hard to beat loop unless you use a MEX solution as "ndfun," which calls the same library BLAS (DGEMM) as mtimes. You have to tradeoff the gains from BLAS versus the overhead of the loop.

Jinhui

[Junior]

Thank you all for the many suggestions.

Jinhui: I took a brief look at multprod. But I could not figure out what is the right way to call the function in this case. I might have missed something with the "IdA" and "IdB". I will take another look.

cheers,
Junior

[Jinhui Bai]

"Junior " <junio...@gmail.com> wrote in message <gq9ro2$fjv$1...@fred.mathworks.com>...

Junior:
In your example, "A = multiprod(B,C)" or "A = multiprod(B,C,[1 2],[1 2])".
Jinhui

[Junior]

Dear All,
I just wanted to report on how I applied the many solutions you suggested for my vectorization problem.
Just to recap, the original problem was to vectorize a the following piece of code:

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

PETER suggested that I use ndfun, which I downloaded but could not run and could not compile
BRUNO suggested a sparse solution for the inner loop only. Does that mean that I have to write an extra loop around the sparse solution?
The other solutions where not restricted to the number of dimensions. MATT suggested that I reverse the order of increments that is instead of going from 1 to n, I go from n to 1 while JINHUI suggest multiprod.
In order to compare the various approaches, I reduced the problem to a one loop problem (unfortunately I have not been able to
implement the mex solution of PETER yet). Because of the instability in the timing results I decided to average several runs

solution loop ::0.0011349
solution loop backwards (MATT) ::0.0010161
solution sparse (BRUNO) ::0.01243
solution multiprod (JINHUI) ::0.014701

Should anyone have any further suggestion or comment, please just shoot.
Thank you for the suggestions.
Junior

The code I ran is below

clear
clc
ra=70;
ca=60;
pa=5;
cb=50;
A=rand(ra,ca,pa);
B=rand(ca,cb,pa);
tloop=0;
tndfun=0;
tmatt=0;
tbruno=0;
tmultprod=0;

N=1000;

for reps=1:N
%% Solution 0
C0=zeros(ra,cb,pa);
tic
for i=1:pa
C0(:,:,i)=A(:,:,i)*B(:,:,i);
end
tloop=tloop+toc;

%% Solution 1: Peter
% tic
% C1 = ndfun('mult', A, B);
% tndfun=tndfun+toc;

%% solution 2: Matt
C2=zeros(ra,cb,pa);
tic
for i=pa:-1:1
C2(:,:,i)=A(:,:,i)*B(:,:,i);
end
tmatt=tmatt+toc;
%% solution 3: Bruno

Ia = repmat(reshape(1:ra*pa,ra,1,pa),[1 ca 1]);
Ja = repmat(reshape(1:ca*pa,1,ca,pa),[ra 1 1]);
Ib = repmat(reshape(1:ca*pa,ca,1,pa),[1 cb 1]);
Jb = repmat(reshape(1:cb*pa,1,cb,pa),[ca 1 1]);
tic
A3 = sparse(Ia(:),Ja(:),A(:));
B3 = sparse(Ib(:),Jb(:),B(:));
C3 = reshape(nonzeros(A3*B3),[ra cb pa]);
tbruno=tbruno+toc;

%% solution 4: Jinhui

tic,C4 = multiprod(A,B);tmultprod=tmultprod+toc;
if reps==1
all([isequal(C4,C0),isequal(C3,C0),isequal(C2,C0)]) % <- Check ,isequal(C1,C0)
end
end

disp(['solution loop ::',num2str(tloop/N)])
% disp(['solution ndfun (PETER) ::',num2str(tndfun/N)])
disp(['solution loop backwards (MATT) ::',num2str(tmatt/N)])
disp(['solution sparse (BRUNO) ::',num2str(tbruno/N)])
disp(['solution multiprod (JINHUI) ::',num2str(tmultprod/N)])

[Bruno Luong]

"Junior " <junio...@gmail.com> wrote in message <gqclqo$c2g$1...@fred.mathworks.com>...

> BRUNO suggested a sparse solution for the inner loop only. Does that mean that I have to write an extra loop around the sparse solution?

No just reshape your 4D in 3D by merging the two last dimensions, then expand it after getting the result.

nb = size(A,3)*size(A,4)
A=reshape(A,[m n nb]);
B=reshape(A,[n p nb]);

... % calculation

C=reshape(C,[m p size(A,3) size(A,4)])

%Bruno

[Bruno Luong]

"Junior " <junio...@gmail.com> wrote in message <gqclqo$c2g$1...@fred.mathworks.com>...

>

> Should anyone have any further suggestion or comment, please just shoot.
> Thank you for the suggestions.
> Junior

My preference is multiprod

I can't compile Peter's ndfun on Vista 64 bit. My lapack library seems to miss two required functions (zemm and demm or something like that).

Bruno

[Jinhui]

On Mar 25, 3:33 am, "Bruno Luong" <b.lu...@fogale.findmycountry>
wrote:
> "Junior " <junior.m...@gmail.com> wrote in message <gqclqo$c2...@fred.mathworks.com>...

>
> > Should anyone have any further suggestion or comment, please just shoot.
> > Thank you for the suggestions.
> > Junior
>
> My preference is multiprod
>
> I can't compile Peter's ndfun on Vista 64 bit. My lapack library seems to miss two required functions (zemm and demm or something like that).
>
> Bruno

Bruno,

It took me two hours to compile ndfun when I first did that two months
ago. I don't know whether this helps you, but this is how I made that
work finally. The key is to link BOTH to LAPACK and BLAS. If you
missed either library, it won't work.

On Windows XP 32 Bit, I did the following with Microsoft Visual C++
2008 Express Edition:
mex ndfun.c -l "C:\Program Files\MATLAB\R2008b\extern\lib
\win32\microsoft\libmwlapack.lib"
-l "C:\Program Files\MATLAB\R2008b\extern\lib\win32\microsoft
\libmwblas.lib"

On Redhat Linux 64 Bit, I did the following with Builtin gcc:
mex ndfun.c -l libmwlapack -l libmwblas

Jinhui

[Bruno Luong]

Thank you Jinhui for the tip of compiling ndfun.

Bruno

[Junior]

Hi Bruno,
Just wanted to warn you that I discovered a possible bug in your formula for computing A(:,:,j)*B(:,:,j). When you reshape the matrices, there is a call to nonzeros, which messes up everything if the product of matrices A and B produces some zeros.

Best,
Junior

[Bruno Luong]

"Junior " <junio...@gmail.com> wrote in message <gqlq3q$a23$1...@fred.mathworks.com>...

Thanks. I got it. It's easy to fix, but not sure it's still fast.

Bruno

[Bruno Luong]

Sparse version, bug fixed:

clear

% Create all data.
m = 3;
n = 4;
p = 2;
nb = 1e5;
A = rand(m,n,nb);
B = rand(n,p,nb);

% Sparse matrix
tic
I = repmat(reshape(1:m,[m 1 1]),[1 n nb]);
J = repmat(reshape(1:n*nb,[1 n nb]),[m 1 1]);
As = sparse(I(:),J(:),A(:),m,n*nb);

I = repmat(reshape(1:n*nb,n,1,nb),[1 p 1]);
J = repmat(reshape(1:p*nb,1,p,nb),[n 1 1]);

Bs = sparse(I(:), J(:), B(:),n*nb,p*nb);
Cs = reshape(full(As*Bs),[m p nb]);
toc %Elapsed time is 0.325316 seconds.

% Bruno

[Matt Fig]

Just a note that this vectorized approach uses a lot of memory, resulting in OUT OF MEMORY errors with fairly small arrays. Also, it is not faster than your For loop in the case of larger arrays. The approach you choose should depend on your expected sizes.

clear all,clc
% Create all data.
m = 10; % Try with m = 20, n = 25! On 32bit winvista, 3GB ram -> OUT OF MEM
n = 15;
p = 10;

nb = 1e5;
A = rand(m,n,nb);
B = rand(n,p,nb);

% Sparse matrix
tic
I = repmat(reshape(1:m,[m 1 1]),[1 n nb]);
J = repmat(reshape(1:n*nb,[1 n nb]),[m 1 1]);
As = sparse(I(:),J(:),A(:),m,n*nb);
I = repmat(reshape(1:n*nb,n,1,nb),[1 p 1]);
J = repmat(reshape(1:p*nb,1,p,nb),[n 1 1]);
Bs = sparse(I(:), J(:), B(:),n*nb,p*nb);
Cs = reshape(full(As*Bs),[m p nb]);

toc % Elapsed time is 4.370445 seconds.

% Bruno

%
tic
C = zeros(m,p,nb);

for k=1:nb
C(:,:,k) = A(:,:,k) * B(:,:,k);
end

toc % Elapsed time is 0.868827 seconds.

>> all(C(:)==Cs(:))
ans =
1

[Junior]

Hi All,

I agree with MATT's point about the memory when the number of pages increases. For the type of application I am interested in, however, it is the size of the matrices that is big. I re-run MATT's example but with different specifications of m,n and nb and my correction
to BRUNO's bug. The general result is that it is still difficult to beat the for loop.

LOOP: Elapsed time is 0.006260 seconds.
INVERTED LOOP (MATT): Elapsed time is 0.004774 seconds.
BRUNO'S SPARSE: Elapsed time is 0.031035 seconds.
MY MODIFICATION TO BRUNO: Elapsed time is 0.010150 seconds.

I'm so sad. I still don't know how to compile ndfun.c. Can anyone bother helping me on this?

Here is the code I ran

clear all,clc
% Create all data.

m = 100;
n = 90;
p = 70;
nb = 4;

A = rand(m,n,nb);
B = rand(n,p,nb);

% Looping

tic
C = zeros(m,p,nb);
for k=1:nb
C(:,:,k) = A(:,:,k) * B(:,:,k);
end
toc

% MATT's original suggestion
tic
Cmatt = zeros(m,p,nb);
for k=nb:-1:1
Cmatt(:,:,k) = A(:,:,k) * B(:,:,k);
end
toc

% BRUNO's correction

tic
I = repmat(reshape(1:m,[m 1 1]),[1 n nb]);
J = repmat(reshape(1:n*nb,[1 n nb]),[m 1 1]);
As = sparse(I(:),J(:),A(:),m,n*nb);
I = repmat(reshape(1:n*nb,n,1,nb),[1 p 1]);
J = repmat(reshape(1:p*nb,1,p,nb),[n 1 1]);
Bs = sparse(I(:), J(:), B(:),n*nb,p*nb);

Cbruno = reshape(full(As*Bs),[m p nb]);
toc

% My suggestion of correction to BRUNO's bug
tic
I=repmat(reshape((1:nb*m),m,nb),n,1);
J=repmat((1:n*nb),m,1);
AA=sparse(I(:),J(:),A(:));
BB=reshape(permute(B,[2,1,3]),p,n*nb);
Cj=permute(reshape(BB*AA',p,m,nb),[2,1,3]);
toc

all(C(:)==Cmatt(:)) && all(C(:)==Cbruno(:)) && all(C(:)==Cj(:))

[Junior]

I finally managed to compile ndfun.c and it wins the day.

Thanks to Peter

[Rodrigo]

How did you get ndfun to work in 64 bits? I tried pointing mex to the blas and lapack libraries included in $MATLABROOT/bin/glnxa64 (I am using the 64bit linux version) and after fixing one or two of the comments it compiles. However, it segfaults whenever I run it.

"Junior " <junio...@gmail.com> wrote in message <grphh1$lne$1...@fred.mathworks.com>...

[Bruno Luong]

"Rodrigo" <guerra.re...@physics.harvard.edu> wrote in message <hnbsg0$rvp$1...@fred.mathworks.com>...

> How did you get ndfun to work in 64 bits? I tried pointing mex to the blas and lapack libraries included in $MATLABROOT/bin/glnxa64 (I am using the 64bit linux version) and after fixing one or two of the comments it compiles. However, it segfaults whenever I run it.

Check out MTIMESX by James Tursa, it's even better than NDFUN for multiple-matrix product.

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

Bruno

[Rodrigo]

I found a rather kludgy way to get it to work by replacing most of the int's with mwSignedIndex's. So in addition to having to point it to -lmwblas and -lmwlapack I also had to add -largeArrayDims. Since I have no idea if this is platfrom independent, if anyone has a better solution please let me know.

I also added an additional routine that mimics cumprod. In other words, given A[i,j,k,...] it produces C[i,j,k,...] st C[:,:,k+1,...]=C[:,:,k,...]*A[:,:,k+1,...]. If anyone is interested I can send it to you or post it.

Rodrigo