The art of vectorizing – Part 2 In a previous post, I introduced you how to the art of vectorize your code. Here I continue in this serie and present some additional tips. I show how to use repmat to process series of elements all at once.

As a neuroscientist, it often came to me that I had to process multiple time traces. Each one of these time traces was an individual trial of the same experiment. So my data matrix Trials was N by T where N is the number of repetition and T the number of time points. A very common thing to do is to substract the average over time from each trial. With for loops you would do :

for i=1:N
AverageValue=mean(Trials(i,:));
Trials(i,:)=Trials(i,:)-AverageValue;
end

To vectorize this code, you first need to get an average matrix. You need to process all the trials at the same time. mean is quite handy for that as mean(Trials,2) will average all the rows and create a column of all the average.
Once you have the average, you face the problem that Trials is a big matrix, so you need to recreate a big average matrix of the same size. repmat is just what you need. It will duplicate your average column over all the T time points.

For instance with A=[1 2 3], repmat(A,2,1) will be

1 2 3
1 2 3

Using repmat, the vectorized version of this for loop is :

AverageValue=mean(Trials,2);
ReplicatedAverage=repmat(AverageValue,1,T);
Trials=Trials-ReplicatedAverage;
This entry was posted in Intermediate, Optimizing your code. Bookmark the permalink.

9 Responses to The art of vectorizing – Part 2

1. Nick S says:

I think you can do better here with bsxfun. At least for this particular example, I get the following:

>> a = rand(1000,10000);
>> tic; mn = mean(a,2); mnrep = repmat(mn, 1,size(a,2)); res1 = a-mnrep; toc
Elapsed time is 0.164605 seconds.
>> tic; res2 = bsxfun(@minus, a, mean(a,2));toc;
Elapsed time is 0.105141 seconds.

(res1 and res2 are identical).

Nick

• Jerome says:

Yes, this is going to be Part 3. I wanted to talk about bsxfun, colfilt, blockproc.
Maybe I’ll spend some time to compare all these functions.

• Jerome says:

Also, I think your measurement is not entirely fair. You should embed this code in an M-file and run it from there.
The JIT is deactivated from the command line, I believe.

2. arijithazra says:

Is there any way we can vectorize matlab for loop involving ode solver at each step and multi-dimensional matrix e.g.

tspan=0:0.01:20;
dw=rand(p,1);
M0=repmat([0 0 1],p,1)’;
for p=1:ns
[t,M(:,:,p)]=ode45(@(t,M) testfun(t,M,dw(p)),tspan,M0(:,p));
end

3. arijithazra says:

where

function dM=testfun(t,M,w1)
M_x=M(1);
M_y=M(2);
M_z=M(3);
dM=[w1*M_y;-w1*M_x+w1*M_z-2*w1*M_y;-w1*M_y-(1-M_z)];

• Jerome says:

I am not sure you can do this straight with matrices. Maybe with Cell if you convert each ‘frame’ to a cell element.
Have you looked into cellfun?

4. John Franklin Ruiz Neira says:

Can some one helpme to vectorizing this section code:

for i=1:alto
for j=1:ancho
E=etiquetas(i,j);
if E~=0
areaZona(E+1)=areaZona(E+1)+1;
end
end
end

5. Math says:

I’m rather new to MatLab and I really need to get used to vectorization. Currently I’m stuck on a problem, where I need to use matrix multiplication in a for loop. I’m wondering if it coud be done without one. The current code:
J=zeros(3,1);
f(size(A,1),size(A,2))=0;
for j=1:size(A,2)
for i=1:size(A,1)
if A(i,j,1)+A(i,j,2)+A(i,j,3)>0
for k=1:7
J(:)=double(A(i,j,:));
aa=(J-M(:,k));
bb=inv(Cov(:,:,k));
ff=(2*pi)^(-3/2)*norm(Cov(:,:,k))^(-1/2)*exp((-1/2)*aa’*bb*aa);
if ff>f(i,j)
f(i,j)=ff;
end
end
end
end
end

A is RGB image (large part of it is black); aa is 3 y 1 vector, bb is 3 by 3 matrix. The current version takes to complete for over 30 min (way too long). Are there a way to use matrix multiplications along one dimension of matrix ‘without’ the loop?

• Jerome says:

Sorry, your posted code is missing the definition of half of the variables. It’s therefore hard to help you.

Is you for loop a bottle neck? Did you time it?

This site uses Akismet to reduce spam. Learn how your comment data is processed.