Efficient way of Vectorization
이전 댓글 표시
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
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.
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
2022년 11월 30일
Mike Croucher
2022년 11월 30일
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
2022년 11월 30일
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
2022년 11월 30일
채택된 답변
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
Nadeem Ahmed
2022년 11월 30일
Mike Croucher
2022년 11월 30일
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
2022년 11월 30일
Bruno Luong
2022년 11월 30일
편집: Bruno Luong
2022년 11월 30일
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
2022년 11월 30일
Bruno Luong
2022년 11월 30일
편집: Bruno Luong
2022년 11월 30일
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
2022년 11월 30일
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
추가 답변 (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
Nadeem Ahmed
2022년 11월 30일
Mike Croucher
2022년 11월 30일
Do it one loop at a time. Depending on the access patter, some loops will be faster...others will not.
Nadeem Ahmed
2022년 11월 30일
Dyuman Joshi
2022년 11월 30일
카테고리
도움말 센터 및 File Exchange에서 Startup and Shutdown에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
