Sunday, May 6, 2007

Getting rid of loops

Programs like Matlab and Octave have problems in doing loops while running user made functions. So it is best to avoid them as much as possible. One simple example is multiplying every element in a matrix A by a corresponding element in matrix B, this could be done by the following loops

for i in [1:size(A(:,1))]
for j in [1:size(A:,2))]
C(i,j) = A(i,j)*B(i,j);
end
end

Now, this would take a long time to run if A and B were really big matrixes (like millions of elements). But if you use the inherent multiplication rules in the programs, something like

C = A.*B;

then this calculation only takes a few seconds.

A friend of mine had a problem the other day where he thought he could not get rid of the for loops and his calculations were taking to long and he didnt even have the patience to see them finished (he was thinking about doing them in C++ or fortran instead) but having a one night look at the problem I found out that his problem could actually be solved a lot faster by using algebraic notation.

His problem was that he had 4 loops, 2 loops to go over each element in a matrix A and two loops to go over a certain amount of neighbors to that element in each of the mentioned loops, something like this

for i in [1:size(A(:,1))]
for j in [1:size(A:,2))]
sa = something1;
sb = something2;
ta = something3;
tb = something4;
for k in [sa:sb]
for l in [ta:tb]
Z = f(i,j,k,l);
end
end
end
end

Now this was a lot more complicated and the function f(i,j,k,l) had a nonlinear factor in it so at first I thought it wasnt possible to simplify this but I was proven wrong fortunately :-)
The trick was that functions can take matrixes as input in Octave and by using the . operator they perform the function on each element in the matrix and not on the matrix as a whole, like is the difference between A*B and A.*B in matrix notation, the first multiplies each column with a row in the other matrix, but the latter multiplies each element together. The same can be done with nonlinear operators like the power operator ^
So for loops like is above where the function is Z = 2^(sqrt(A(i,k)^2 + A(j,l)^2))) can be simplified by using the meshgrid command

x = [sa:sb]
y = [ta:tb]
[X,Y] = meshgrid(x,y)
Z = 2^(sqrt(X.^2 + Y.^2))

does the same thing except much faster. I didnt manage to get rid of the first two loops but he says that his runtime went down from infinite to about 5 minutes, thats pretty impressive I think.

Hope this can help someone who is learning Matlab or Octave

No comments: