why vectorization becomes slower?

I was trying to compute a simple function by 'for' loop but when I vectorized the same function as shown below in 'VECTORISED FORM' instead of improvement in CPU time it becomes slower than 'for' loop, I don't exactly know the reason for this because I am new in vectorization therefore don't even know whether I did vectorization right or wrong

SIMPLE FORM

for i=2:n+1
    for j=2:m+1
        f1(i,j)=  -((+ rho*(((((uyt(i,j) - uyt(i-1,j))/(2*a))+((vxt(i,j) - vxt(i,j-1))/(2*a))))/(2*tao))) + ...
            (rho*(((uytold(i,j)-uytold(i-1,j))/(2*a))*((uyt(i,j)-uyt(i-1,j))/(2*a)))) + ...             
            (rho*(((vyt(i,j)-vyt(i-1,j))/(2*a))*((uxtold(i,j)-uxtold(i,j-1))/(2*b)))) + (rho*(((vytold(i,j)-vytold(i-1,j))/(2*a))*((uxt(i,j)-uxt(i,j-1))/(2*b)))) + ...
            (rho*(((vxtold(i,j)-vxtold(i,j-1))/(2*b))*((vxt(i,j)-vxt(i,j-1))/(2*b)))));                  
    end
end

VECTORISED FORM

f1(2:n+1,2:m+1)=  -((+ rho.*(((((uyt(2:n+1,2:m+1) - uyt(1:n,2:m+1))./(2.*a))+((vxt(2:n+1,2:m+1) - vxt(2:n+1,1:m))./(2.*a))))./(2.*tao))) + ... 
 (rho.*(((uytold(2:n+1,2:m+1)-uytold(1:n,2:m+1))./(2.*a)).*((uyt(2:n+1,2:m+1)-uyt(1:n,2:m+1))./(2.*a)))) + ...              
 (rho.*(((vyt(2:n+1,2:m+1)-vyt(1:n,2:m+1))./(2.*a)).*((uxtold(2:n+1,2:m+1)-uxtold(2:n+1,1:m))./(2.*b)))) + (rho.*(((vytold(2:n+1,2:m+1)-vytold(1:n,2:m+1))./(2.*a)).*((uxt(2:n+1,2:m+1)-uxt(2:n+1,1:m))./(2.*b)))) + ...
(rho.*(((vxtold(2:n+1,2:m+1)-vxtold(2:n+1,1:m))./(2.*b)).*((vxt(2:n+1,2:m+1)-vxt(2:n+1,1:m))./(2.*b)))));      

Answer 1

[Matt J]

One reason is that you are extracting the same submatrices multiple times. For example, the expression uyt(2:n+1,2:m+1) appears multiple times in your vectorized code. This operation allocates memory for a new matrix every time you do it. You should really go through your expression and pull out repeated instances of sub-expressions and precompute them.

Also, it looks like you can replace certain expressions like uyt(2:n+1,2:m+1) - uyt(1:n,2:m+1) with diff(). If uyt is of size (n+1)x(m+1), then the latter will be faster, e.g.,

m=2000;
n=2000;

uyt=rand(n+1,m+1);

tic
d_uyt = uyt(2:n+1,2:m+1) - uyt(1:n,2:m+1);
toc

Elapsed time is 0.053775 seconds.

tic
d_uyt = diff(uyt,1,1);
toc

Elapsed time is 0.027445 seconds.

[Nadeem Ahmed]

Thank you Sir, I can try this 'diff' operation but I can not do anything with the expression it should be like that because it is a non linear function and I dont understand the 'Precompute' thing can you please give an example because at the end I have to put all the things in the same expression.

[Matt J]

I guess you figured it out, since you accepted my answer?

[Nadeem Ahmed]

We can also see the simpler case

n=100;
m=100;
F53=rand(n+1,m+1);
F54=rand(n+1,m+1);
F57=rand(n+1,m+1);
F74=rand(n+1,m+1);
F75=rand(n+1,m+1);
F77=rand(n+1,m+1);
d1=zeros(1,n-1);
d2=zeros(1,n-1);

%NORMAL FORM
tic
for j=2:m+1
    for i=2:n
        d2(i-1)=-F57(i,j) + (F53(i,j)*F74(i,j))/F77(i,j) + (F54(i,j)*F75(i+1,j))/F77(i+1,j);
    end
end
toc  %Elapsed time is 0.000938 seconds.

%VECTORIZED FORM
tic
for j=2:m+1
    d1(1:n-1)=-F57(2:n,j) + (F53(2:n,j).*F74(2:n,j))./F77(2:n,j) + (F54(2:n,j).*F75(3:n+1,j))./F77(3:n+1,j);
end
toc  %Elapsed time is 0.001801 seconds.

Why Vectorised form takes more time and what I did wrong in this simpler case?

[Matt J]

Again, you have lots of non-scalar sub-matrix extractions like F57(2:n,j). These are costly. Look at the following slightly different comparison.

n=2000;
m=2000;
F53=rand(n+1,m+1);
F54=rand(n+1,m+1);
F57=rand(n+1,m+1);
F74=rand(n+1,m+1);
F75=rand(n+1,m+1);
F77=rand(n+1,m+1);
d1=zeros(1,n-1);
d2=zeros(1,n-1);

%NORMAL FORM
tic
for j=2:m+1
    for i=2:n
        d2(i-1)=-F57(i,j) + (F53(i,j)*F74(i,j))/F77(i,j) + (F54(i,j)*F75(i+1,j))/F77(i+1,j);
    end
end

toc 

Elapsed time is 0.170808 seconds.

%VECTORIZED FORM
F53=F53(2:n,:);
F54=F54(2:n,:);
F57=F57(2:n,:);
F74=F74(2:n,:);
F75=F75(3:n+1,:);
tmp1=F77(2:n,:);
tmp2=F77(3:n+1,:);

tic
d1=-F57 + (F53.*F74)./tmp1 + (F54.*F75)./tmp2;
toc 

Elapsed time is 0.019281 seconds.

[Nadeem Ahmed]

I totally get it but the main problem is that you define all the submatrices outside and doesn't count the time for them thats why it is showing less time. In my case I will call the function (vectorisation is inside function) in my main program and the cost of whole main program will include these submatricies cost also therefore we can not neglect them.

[Nadeem Ahmed]

n=2000;
m=2000;
F53=rand(n+1,m+1);
F54=rand(n+1,m+1);
F57=rand(n+1,m+1);
F74=rand(n+1,m+1);
F75=rand(n+1,m+1);
F77=rand(n+1,m+1);
d1=zeros(1,n-1);
d2=zeros(1,n-1);

%VECTORIZED FORM
tic
F53=F53(2:n,:);
F54=F54(2:n,:);
F57=F57(2:n,:);
F74=F74(2:n,:);
F75=F75(3:n+1,:);
tmp1=F77(2:n,:);
tmp2=F77(3:n+1,:);

d1=-F57 + (F53.*F74)./tmp1 + (F54.*F75)./tmp2;
toc 

Elapsed time is 0.191028 seconds.

[Matt J]

Yes, hopefully you now see my point. To gain the benefits of vectorization, you must avoid the sub-matrix extraction operations. They are costly.

[Matt J]

So, you have a few options.

(1) Accept the performance of the original loop. Are you sure it is the bottleneck in your code?

(2) Implement the loop in a MEX file. In a MEX, you are not forced to extract data copies of sub-matrices. This can be worthwhile if the code is truly a critical bottleneck (see also option 1).

(3) When you create your matrices, try to avoid embedding the data blocks you'll need later inside larger matrices. In the code you've shown, for example, the first and last row of F53 are never used. Do they need to be there? If you need them, could they be kept in their own separate variables?

[Nadeem Ahmed]

Thank you for your suggestions.