The art of vectorizing – Part 1

Vectors and matrices are at the heart of any Matlab program. Matlab is extremely efficient at doing any operation with these. Vectorizing is the art of transforming a calculation done element by element into an operation on vectors. Here I start a serie of posts on how to vectorize your code.

I had a very good teacher in chemistry when I was in college. He used to say something like this :

To understand chemistry is easy, you just have to understand the very few building blocks on which all this science is based upon. And they don’t take more space than a metro ticket…

Programming is similar, I believe. They are some important principles and all the rest is just addon knowledge built on these principles. If you were to ask me, vectorizing is one of the few Matlab core principles that should sit on this metro ticket. Some day I will make a post to actually draw this metro ticket in detail…

I already talked about why you should vectorize in numerous occasions here, here and here. So today I am going to dive straight into the hard core. How to vectorize?

They are many ways to process data in blocks and I want to cover all the tips and tricks I have used in my Matlab life so I will make as many posts as needed on this.

But I guess to start, we should talk about what I consider the most important Matlab operator : the point.

Such a small thing and yet I remember my first days with Matlab. This point was not intuitive. Part of the reason is because I was not taught Mathematic this way.

Let’s assume you have two matrices A and B of the same size. There are two ways to multiply them element by element :

A=rand(100,1);
B=rand(100,1);
for i=1:100
C(i)=A(i)*B(i);
end

Or you could use the magical . operator. The point actually modify the multiplication operator so that it operates element by element, instead of doing a complicated matrix multiplication.

A=rand(100,1);
B=rand(100,1);
C=A.*B;

Please note that the small point. A.*B is not equal to A*B. Here A*B will give you an error because you are not using properly the matrix multiplication (mathematically correct operations would be A’*B or A*B’ here where ‘ is the transpose). So the point MODIFIES the definition of * to make a multiplication of the entire A and B, element by element.

If you were to ask me, I believe this should be the default behavior of matrix multiplication but I can’t restart Mathematic history…

The point is (I couldn’t resist…) that this is EXTREMELY usefull and is probably one of the simplest form of vectorization you could do. The vectorize version A.*B is much faster than the for loop as Matlab operates on the entire matrix at once via its built-in routines.

16 Responses to The art of vectorizing – Part 1

Hi, I would like you to know that I have found your website really useful, and even enjoyable. over the years, I look topics up via Google and I’ve come across your website so many times that I recognize it without reading anything on the page! I’m finally going to bookmark your website, but, I’ve been learning from you all along. I’m an undergrad by the way; Applied Math and Materials Engineering major – thank you for your hard work.

I have always heard vectorizing is the way go but never tried too hard to go out of my way to make sure I was doing it — I guess because what I was doing just wasn’t that time sensitive. Today I thought I’d try out your suggestion. However I’m seeing some unusual results.

Can you think of any reason why if I use tic…toc around the 2 examples you give above, the for loop comes out on top?

Here’s the code:

A=rand(100,1);
B=rand(100,1);

tic
for i=1:100
if B(i)>0.5
C(i)=A(i)^2;
else
C(i)=exp(B(i));
end
end
toc

I think you are seeing a mixture of preallocation time + JIT (for Just In Time) compilation. JIT can speed up computation of for loop quite significantly. The JIT is especially efficient on these simple for loops.

Thanks for the reply. That does make a lot of sense. I think I did run it a couple times before the numbers in my original post and I had about the JIT improving things. I have an older version of MATLAB kicking around on another computer that I may try this on out of curiosity, as well.

Thanks again for the reply and for this site (I just stumbled upon it and plan on digging around a bit more).
Cheers

I tried the following code, but i ran into trouble. A as a vector of real numbers, some of which are NaN, because sometimes a division by zero occurs when calculating A. I want to replace these NaN’s with zeros.

NANs=(isnan(A));
A=(~NANs).*A + (NANs).*0;

It doesn’t seem to work, because, 0 * NaN = NaN.

Do you have a clever way to do this in a vectorized manner, or do I have to resort to a for-loop?

Aah! I didn’t expect that to work, because I know e.g. “A([0 1])=0” returns an error because the first index is zero. However your suggestion works because isnan() returns an array of logicals the same size as A, and not an array of numbers.

I don’t think you should vectorize every loops. First, sometimes, it is not possible or easy to do so. Second, sometimes it goes against clarity and provide very little improvement in efficiency. Use common sense and proper evaluation of your bottle necks.

my GUESS is, it is possible if all variables are either matrices of the same size as “img”, or scalars. So if “noupdateMask” is the same size as as “img(:,:,i)” (i.e., two-dimensional), then you would have to make it three dimensional with repmat(). This way you avoid ambiguities where MATLAB doesn’t know how to apply the maths to the different dimensions.

Do note however, that you exchange CPU time for memory usage. Every intermediate “mask”, as well as “noupdateMask”, and “refLoc” will become size(img,3) times larger in memory. Which might slow things down too for large data sets (such as a movie?), or even freeze up the PC. (Experience!)

You already have this vectorized along two dimensions, which is already pretty good, I’d say.

Thank you very much Pieter van Vugt! I find your response good and i thought the same thing!!
For information:“mask”, as well as “noupdateMask”, and “refLoc” have the same dimension as «img»!!
I am waiting for Jerome’s confirmation 🙂 and it will be pleasant 🙂

Hey Jerome. Just a question about your last example. It seems counter intuitive to do it the way you do this, since the for loop would take advantage of temporal memory? (something something more cache hits) Am I wrong in my line of reasoning?