Efficient way of Vectorization
Mostra commenti meno recenti
Hello, I searched everywhere for the efficient explaination of vectorization, I would like to know how can we use the technique of vectorization efficiently if we have this kind of problem
clc
clear
close all
n=1000;
C1=zeros(n,n);
C2=zeros(n,n);
A=rand(n,n);
B=rand(n,n);
tic
for i=2:n-1
for j=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
toc
%VECTORIZATION
tic
C2(2:n-1,2:n-1)=(A(2:n-1,2:n-1).*B(2:n-1,1:n-2) + A(1:n-2,2:n-1).*B(3:n,1:n-2))./(A(2:n-1,3:n).*B(3:n,2:n-1));
toc;
norm(C1-C2)
This is a very basic example, although it is showing the improvement after vectorization but not that enough. If I make more divison and multiplication in the same function, "vectorization" will become even worse than "for loop". If anybody have any suggestion regarding this, it would be very helpful for me.
8 Commenti
You can initialize the index vector, so you don't have to define it everytime.
clc
clear
close all
n=1000;
C1=zeros(n,n);
C2=zeros(n,n);
A=rand(n,n);
B=rand(n,n);
tic
for i=2:n-1
for j=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
toc
%VECTORIZATION
tic
C2(2:n-1,2:n-1)=(A(2:n-1,2:n-1).*B(2:n-1,1:n-2) + A(1:n-2,2:n-1).*B(3:n,1:n-2))./(A(2:n-1,3:n).*B(3:n,2:n-1));
toc
%Vectorization with index vector
v=2:n-1;
tic
C3(v,v)=(A(v,v).*B(v,v-1) + A(v-1,v).*B(v+1,v-1))./(A(v,v+1).*B(v+1,v));
toc
Though, from what I know, using tic-toc is not the most accurate method, as the run-time depends on many other factors. You can either run the same code repeatedly using a for loop and get the mean time or use timeit() (which outputs median time)
However, for crude comparison, you will find that vectorization with index vector is faster than plain vectorization by a goood margin.
Nadeem Ahmed
il 30 Nov 2022
Modificato: Nadeem Ahmed
il 30 Nov 2022
As I mentioned earlier, results from tic-toc can vary
n=10000;
C1=zeros(n,n);
C2=zeros(n,n);
C3=zeros(n,n);
A=rand(n,n);
B=rand(n,n);
tic
for i=2:n-1
for j=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
toc
%VECTORIZATION
tic
C2(2:n-1,2:n-1)=(A(2:n-1,2:n-1).*B(2:n-1,1:n-2) + A(1:n-2,2:n-1).*B(3:n,1:n-2))./(A(2:n-1,3:n).*B(3:n,2:n-1));
toc
%Vectorization with index vector
v=2:n-1;
tic
C3(v,v)=(A(v,v).*B(v,v-1) + A(v-1,v).*B(v+1,v-1))./(A(v,v+1).*B(v+1,v));
toc
However, without knowing what do you exactly want to do, it's difficult to say how you can vectorize efficiently further. And with regard to this particular example, it seems difficult to increase efficiency.
Nadeem Ahmed
il 30 Nov 2022
Mike Croucher
il 30 Nov 2022
Can you give an example of how to call this and how long it takes for you? I want to run
tic
MAIN(v1,v2,v3,v4,v5,v6,v7,v8,l11,l22,l33,l44,d11,d22,d33,d44,u11,u22,u33,u44,F11,F13,F15,F16,F17,F21,F23,F25,F26,F27,F33p,F34p,F37p,F53p,F54p,F57p,F71p,F72p,F74p,F75p,F77p,uxtold,uytold,vxtold,vytold,Pxtold,Pytold,tao,a,b,m,n,P_inv)
toc
but I don't know what v1,v2,v3 etc are for an example you care about.
Mike Croucher
il 30 Nov 2022
Thanks. So for N,M=50, the code runs in 0.01 seconds on my machine.
Increasing to N,M=100, the code runs in 0.22 seconds
Trying N,M=200, I run out of memory on my 32Gb laptop.
What values of N and M are you interested in and how fast do you need the code to be?
Nadeem Ahmed
il 30 Nov 2022
Risposta accettata
Unfortunately, this is a situation where the for loop is the fastest option. This is because vectorized solution does much more memory allocation than it should. I have raised this issue with MathWorks staff, but am not sure what is being done on it.
function test
n=1000;
C1=zeros(n,n);
C2=zeros(n,n);
A=rand(n,n);
B=rand(n,n);
timeit(@()method1)
timeit(@()method2)
ans =
0.0161
ans =
0.0210
function method1
for i=2:n-1
for j=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
end
function method2
C2(2:n-1,2:n-1)=(A(2:n-1,2:n-1).*B(2:n-1,1:n-2) + A(1:n-2,2:n-1).*B(3:n,1:n-2))./(A(2:n-1,3:n).*B(3:n,2:n-1));
end
end
9 Commenti
Nadeem Ahmed
il 30 Nov 2022
Mike Croucher
il 30 Nov 2022
MATLAB now does JIT (Just In Time) compilation and its pretty good. This is one reason why vectorisation isn't always faster...although vectorisation can take advantage of multithreading if the vectors are large enough.
If you really want to try C++, you could try MATLAB coder to do it automatically.
Nadeem Ahmed
il 30 Nov 2022
Bruno Luong
il 30 Nov 2022
Modificato: Bruno Luong
il 30 Nov 2022
MATLAB now does JIT (Just In Time) compilation and its pretty good. This is one reason why vectorisation isn't always faster
But it is puzzling that the JIT does not also perform well enough to optimize the vectorized version, so that it is competitive with the for-loop. In this example, the vectorized version is wasting time needlessly because subsref operations on A and B on the right hand side of,
C2(2:n-1,2:n-1)=(A(2:n-1,2:n-1).*B(2:n-1,1:n-2) + A(1:n-2,2:n-1).*B(3:n,1:n-2))./(A(2:n-1,3:n).*B(3:n,2:n-1));
are allocating memory for copies of the sub-matrices like A(2:n-1,2:n-1). This is unnecesary because all operations on A and B here are read-only, and one wonders why the JIT isn't able to recognize this.
Mike Croucher
il 30 Nov 2022
Bruno Luong
il 30 Nov 2022
Modificato: Bruno Luong
il 30 Nov 2022
I don't think the problem is allocating memory, but actually indexing with truncation index, which requires elements in memory to be rearranged.
I'm not surprised that to make a vectorize code as fast as the for-loop requires a big development of the internal engine (for instant using meta data that describe subarray of an array without copying the data).
Indexing is always the bottleneck of MATLAB.
Matt J
il 30 Nov 2022
I don't think the problem is allocating memory, but actually indexing with truncation index
Not sure what a "truncation index" refers to here. In any case, the subsref operations are definitely to blame, since when we revise the test with the indexing done offline, the vectorized version is much more competitive with the loops:
function test
n=1000;
C1=zeros(n,n);
C2=zeros(n,n);
A=rand(n,n);
B=rand(n,n);
[Q1,Q2,Q3,Q4,Q5,Q6]=...
deal( A(2:n-1,2:n-1) , B(2:n-1,1:n-2), A(1:n-2,2:n-1),...
B(3:n,1:n-2), A(2:n-1,3:n), B(3:n,2:n-1) );
timeit(@()method1)
timeit(@()method2)
ans =
0.0149
ans =
0.0051
function method1
for i=2:n-1
for j=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
end
function method2
C2(2:n-1,2:n-1)=(Q1.*Q2 + Q3.*Q4)./(Q5.*Q6);
end
end
Più risposte (1)
Switch the order of the loops around. It will be faster because you'll be operating on the matrix column-wise
test
function test
n=2000;
C1=zeros(n,n);
C2=zeros(n,n);
A=rand(n,n);
B=rand(n,n);
disp('loops')
timeit(@()loops)
disp('loops 2')
timeit(@()loops2)
disp('vector')
timeit(@()vector)
function loops
for i=2:n-1
for j=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
end
function loops2
for j=2:n-1
for i=2:n-1
C1(i,j) = (A(i,j)*B(i,j-1) + A(i-1,j)*B(i+1,j-1))/(A(i,j+1)*B(i+1,j));
end
end
end
function vector
C2(2:n-1,2:n-1)=(A(2:n-1,2:n-1).*B(2:n-1,1:n-2) + A(1:n-2,2:n-1).*B(3:n,1:n-2))./(A(2:n-1,3:n).*B(3:n,2:n-1));
end
end
4 Commenti
Nadeem Ahmed
il 30 Nov 2022
Mike Croucher
il 30 Nov 2022
Do it one loop at a time. Depending on the access patter, some loops will be faster...others will not.
Nadeem Ahmed
il 30 Nov 2022
Dyuman Joshi
il 30 Nov 2022
Categorie
Scopri di più su Startup and Shutdown in Centro assistenza e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
