The art of vectorizing – Part 3

In this post, I keep going in my series of posts on vectorizing. Following the last post on using repmat to avoid unnecessary for loops, I introduce bsxfun which is a faster alternative to repmat.

As mentioned recently by Nick, repmat is a nice way to duplicate some arrays around multiple dimensions and avoid some for loops. It is quite a flexible and general approach  that not only limits to processing row by row (or column  by column). You can also replicate various subarrays. However, if your particular application is to process an array either element by element or along one of its dimension, then bsxfun is a faster alternative, as pointed by both Nick on this blog and Loren on her own, especially since bsxfun is now multithreaded in Matlab. Another advantage of bsxfun is that it does not require to allocate an intermediate large matrix. If your particular application is demanding on memory, than it is definitely going to be faster for you.

Besides, there are really no reason not to use it as it is an easier alternative in most cases. Our previous example was, given a NxT matrix Trials :

N=100;
T=1000;
Trials=rand(N,T);
for i=1:N
   AverageValue=mean(Trials(i,:));
   Trials(i,:)=Trials(i,:)-AverageValue;
end

The equivalent with bsxfun is :

N=100;
T=1000;
Trials=rand(N,T);
AverageValue=mean(Trials,2);
Trials=bsxfun(@minus,Trials,AverageValue);

As indicated on the doc, bsxfun applies @minus to 2 matrices A and B (here A is Trials and B is AverageValue). All dimensions of A and B must be the same or singletons. bsxfun will automatically duplicates arrays (in a way automatically applies repmat) so that the operation (here minus) can be applied on all single elements. ‘@’ is used to give a function handle. A function handle is like a pointer that gives to Matlab the address to a particular function that can take both A and B as inputs. There are a number of built-in handles that you can pass to bsxfun. Most standard operations are there. You can also create your own as long as you respect its syntax.

Related Posts

This entry was posted in Intermediate, Optimizing your code. Bookmark the permalink.

6 Responses to The art of vectorizing – Part 3

  1. panos says:

    Hi, I’d like to ask wether it’s possible to vectorize a for loop that looks like this:

    x=[1:100];
    a=3;
    b=5;
    y=zeros(0,length(x));

    for n=1:length(x)
    y(n)=x(n)+a*y(n-b)-a*x(n-b) % or something similar to this
    end

    thanks in advance

    • Jerome says:

      Like this, I don’t see any very obvious way to vectorize that, anybody?. I would have say use the ‘filter’ function but your calls to x(n) make this not possible.
      Is this really slow?
      Like that, if n is only 100, it should be fine.
      Sometimes, you should trust the JIT (Check my post on Matlab and for loops) and keep coding in the more natural way. It is sometimes easier to read and modify.

  2. Jay K. says:

    My impression is that the original poster (panos) is testing you, Jerome. I think he may already know an answer to his question. This is just my guess, though. Vectorization is not always obvious, because the routine to be vectorized often requires a mathematical description for better understanding. Otherwise, it is very difficult to see what the code is trying to achieve. The above is such an example. Once understood mathematically, it can be easily vectorized. Since I don’t have time to fully describe the purpose of the code, I will simply give an answer.

    The vectorized code is:

    yy = x + reshape(exp(bsxfun(@plus, log(y(1:b)-x(1:b))’, log(a)*(0:(length(x)/b-1)))), 1,length(x));

    Note. y(1:b) is should be given in advance. Otherwise, the routine does not make sense. If you run the original code, it will give an error message. You need to define y(1:b) first, and loop over (b+1):length(x). Then, my one-line code will be equivalent to that.

    • Jerome says:

      Thanks Jay for posting an answer.

      If the point was to test my expertise: I don’t intend to be the world expert on vectorizing. Just trying to help folks to get started with Matlab in every way I can. Teaching is the best way to learn and is a very nice social endeavor that I enjoy. So far this blog seemed to have been useful as I can tell from the traffic it generates. I will try to keep the quality as high as I can and highly value all of your feedbacks to do so.

      • Jay K. says:

        Hi Jerome,
        Yes, I like your blog. I enjoyed reading your article about JIT and learned new things. Matlab keeps improving. Mathworks is quick to adopt new technology. For example, I was pleasantly surprised to learn that some new parallel random number generators, proposed in 2011 academic papers, were implemented in Matlab 2012b. Do you know anything about their plan to add multiple-GPU functionality (i.e., computing with all availalbe GPU cards at the same time)?

        • Jerome says:

          Thanks.
          I did play with single GPU processing card in the past. Not extensively yet. I plan to do so in the future. I don’t know about multi-GPU capabilities of Matlab. Maybe Mathworks folks will be able to answer you on this. So far, my time bottleneck has been code design and writing, not running time so I have been focusing on improving that.

Leave a Reply