Question on speed of for loop and repmat/repelem 'vectorization'

I am studying a scheduling problem which concerns on assigning series of jobs to a number of operators with multiple time windows. Suppose the assignment and time consumption of each job is known and determined. We need find out the earliest time window of each operator that is long enough to schedule the task assigned to it. Let NJ and NTW be the total number of jobs (operators) and number of time windows of each operator.

I wrote the following codes to test the performance:

NJ=100;
NTW=100;
Task=[(1:NJ)',4+randi(5,NJ,1)];
TW=[repelem((1:NJ)',NTW,1),10*repmat((1:NJ)'-1,NTW,1),10*repmat((1:NJ)'-1,NTW,1)+5+randi(4,NJ*NTW,1)];
TW(NTW:NTW:NJ*NTW,3)=TW(NTW:NTW:NJ*NTW,2)+10;

tic
A=[];
for j=1:NJ
    TWj=TW(TW(:,1)==j,:);
    stInd=find(TWj(:,3)-TWj(:,2)>=Task(j,2),1);
    st=TWj(stInd,2);
    A=[A;j,st,st+Task(j,2)];
end
toc

I attempted to vectorize the codes using repmat and repelem such that the for-loop can be avoided:

tic
M=[repelem(Task,NTW*NJ,1),repmat(TW,NJ,1)];
M(M(:,1)~=M(:,3),:)=[];
M(M(:,5)-M(:,4)<M(:,2),:)=[];
st=accumarray(M(:,1),M(:,4),[],@min);
A=[Task(:,1),st+Task(:,2)];
toc

However the performance is in fact worse than using for-loop. I searched on Google and found this website: http://www.vincentcheung.ca/research/matlabindexrepmat.html

According to it, repmat is slow as it is actually a Matlab script with possibly large overhead. The author also provided a indexing-based method to replicate the matrix without using repmat. I modified the codes and wrote a function for element-wise replication, which does the same thing as repelem. The codes are as follows:

tic
M=[relm(Task,NTW*NJ,1),rmat(TW,NJ,1)];
M(M(:,1)~=M(:,3),:)=[];
M(M(:,5)-M(:,4)<M(:,2),:)=[];
st=accumarray(M(:,1),M(:,4),[],@min);
A=[Task(:,1),st+Task(:,2)];
toc

function B=rmat(A,i,j)
B=A((1:size(A,1))'*ones(1,i),(1:size(A,2))'*ones(1,j));
end

function B=relm(A,i,j)
rowInd=(1:size(A,1))'*ones(1,i);
rowInd=reshape(rowInd',numel(rowInd),1);
colInd=(1:size(A,2))'*ones(1,j);
colInd=reshape(colInd',numel(colInd),1);
B=A(rowInd,colInd);
end

My finding is that the time spent of these functions is even worse that the built-in repmat and repelem. The running results on my computer is as follows. I am using Matlab R2018a:

>> test
Elapsed time is 0.039591 seconds.
Elapsed time is 0.670042 seconds.
Elapsed time is 1.038526 seconds.

The difference is even more significant with larger problem scale.

Now I wonder that is it still wise to 'vectorize' codes instead of using for loops in current version of Matlab? The performance of loops has become better and better through several recent versions.

[KALYAN ACHARJYA]

I am also waiting for experts comments.

Answer 1

[OCDER]

I don't get what you are trying to compare, as the output "A" is different. Placements of tic/toc are also unfair as accumarray is included in one, but excluded in the other.

Here's a cleaner example of difference in vectorized vs for loops.

N = rand(2000);

tic
M1 = zeros(size(N));
for j = 1:numel(N)
    M1(j) = N(j)^2; %Not vectorized
end
toc  %Elapsed time is 0.043167 seconds.

tic
M2 = N.*N; %Vectorized
toc  %Elapsed time is 0.014563 seconds.

isequal(M1, M2); %make sure it returns 1, otherwise not a FAIR comparison

%Also, consider using the timeit function instead of tic/toc for comparison of function speeds

Here's the explanation why the two examples are slow.

tic
A=[]; %Do NOT grow A with the loop. Preallocate A for increased speed.
for j=1:NJ
    TWj=TW(TW(:,1)==j,:);
    stInd=find(TWj(:,3)-TWj(:,2)>=Task(j,2),1);
    st=TWj(stInd,2);
    A=[A;j,st,st+Task(j,2)]; %This is bad. Causes delete + copy of A to grow it.
end
toc

tic
M=[repelem(Task,NTW*NJ,1),repmat(TW,NJ,1)];  %This is creating a fairly large matrix...
%M(M(:,1)~=M(:,3),:)=[]; %Deleting elements is slower than extracting them.
%M(M(:,5)-M(:,4)<M(:,2),:)=[]; %Same issue as above. Also add parentheses for clarity in index. 
M = M(M(:,1)==M(:,3),:);
M = M((M(:,5)-M(:,4))>=M(:,2),:);
st=accumarray(M(:,1),M(:,4),[],@min); %You included the accumarry for this tic/toc, but not the other case...
A=[Task(:,1),st+Task(:,2)];
toc

For speed issues, begin by reading articles about coding performance and create FAIR comparisons of vectorized vs non-vectorized codes.

https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

https://www.mathworks.com/company/newsletters/articles/programming-patterns-maximizing-code-performance-by-optimizing-memory-access.html

Answer 2

[Bruno Luong]

Please Try this (I hope there is no mistake):

NJ=100;
NTW=100;
Task=[(1:NJ)',4+randi(5,NJ,1)];
TW=[repelem((1:NJ)',NTW,1),...
   10*repmat((1:NJ)'-1,NTW,1),...
   10*repmat((1:NJ)'-1,NTW,1)+5+randi(4,NJ*NTW,1)];
TW(NTW:NTW:NJ*NTW,3)=TW(NTW:NTW:NJ*NTW,2)+10;

tic
A=[];
for j=1:NJ
    TWj=TW(TW(:,1)==j,:);
    stInd=find(TWj(:,3)-TWj(:,2)>=Task(j,2),1);
    st=TWj(stInd,2);
    A=[A;j,st,st+Task(j,2)];
end
toc % Elapsed time is 0.009305 seconds.

% Vectorize way
tic
TWD = [TW(:,1), TW(:,3)-TW(:,2)];
[~,is] = sortrows([Task;TWD]);
b = is > NJ;
i = find(b);
j = find(~b);
[~,k] = histc(j,[-Inf;i]);
i = is(i(k))-NJ;
j = is(j);
st = TW(i,2);
A = [j,st,st+Task(j,2)];
toc % Elapsed time is 0.002887 seconds.
% check the result
all(TW(i,3)-TW(i,2)>=Task(j,2) & TW(i,1)==j)