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;

*Related*

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

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.

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.

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

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)];

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?

Can some one helpme to vectorizing this section code:

areaZona=zeros(cantidadEtiquetas+1,1);

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

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?

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?

Hi Jerome,

Thanks for your articles on vectorization. The Part 1 was particularly helpful to me. However, I am struggling with vectorizing this parfor loop. I want to completely remove the parfor loop from the code as it is taking a long time to execute when n is large. Please see the code pasted below. I will appreciate any tips/advice/help you or any other person in this forum can give me on this. Many thanks in advance.

% Initialization and precomputations

% w is an n x 1 vector

% beta: any number larger than 0. Usually set to 1.

lambda = zeros(n,1);

x = w;

y = w;

v = lambda – (rho.*y);

rhow = rho.*w;

n = length(w);

parfor i = 1 : n

if w(i) >= 0

if v(i) < -rhow(i) – beta – 1

x(i) = (-beta -1 -v(i))./rho;

elseif (-rhow(i) – beta – 1 <= v(i)) && (v(i) <= -rhow(i) + beta – 1)

x(i) = w(i);

elseif (-rhow(i) + beta – 1 < v(i)) && (v(i) < beta – 1)

x(i) = (beta – 1 -v(i))./rho;

elseif (beta – 1 <= v(i)) && (v(i) <= beta + 1)

x(i) = 0;

else

x(i) = (beta + 1 – v(i))./rho;

end

else

if v(i) < -beta -1

x(i) = (-beta -1 – v(i))./rho;

elseif (-beta – 1 <= v(i) )&& (v(i) <= -beta + 1)

x(i) = 0;

elseif (-beta + 1 < v(i)) && (v(i) < -rhow(i) – beta + 1)

x(i) = (-beta + 1 – v(i))./rho;

elseif (-rhow(i) – beta + 1 <= v(i)) && (v(i) <= -rhow(i) + beta + 1)

x(i) = w(i);

else

x(i) = (beta + 1 – v(i))./rho;

end

end

end