2

I am writing a Matlab application that computes some sort of matrix. I'm trying to replace the for loops in my program with vector-based calculations, however, I'm stuck.

So far, I have figured that this simple sample part of code: 

kernel = zeros(5,5,5);
offset = 3;
sig=1;

for x=-2:2
    for y=-2:2
        for z=-2:2
            IN = x.^2/(sig^2) + y.^2/(sig^2) + z.^2/(sig^2); 
            kernel(offset+x, offset+y, offset+z) = exp(-IN/2);
        end
    end
end

can be replaced with such a construction:

[x,y,z] = ndgrid(-2:2,-2:2,-2:2);
IN = x.^2/(sig^2) + y.^2/(sig^2) + z.^2/(sig^2); 
kernel = exp(-IN/2);

and give the same results. But what if I need to make some small alteration:

kernel = zeros(5,5,5);
offset = 3;   
sig=1;

%sample 3x3 matrix
R=magic(3)/10; 

for x=-2:2
    for y=-2:2
        for z=-2:2

            % calculation of new x, y, z
            point = [x,y,z]*R;   

            IN = (point(1)^2 )/(sig^2) + (point(2)^2)/(sig^2) + (point(3)^2)/(sig^2);
            kernel(offset+x, offset+y, offset+z) = exp(-IN/2);
        end
    end
end

How can I speed up this construction? Can it be easily vectorized? I'm quite new to Matlab, so I would appreciate any help. Thanks a lot!

user3299285
  • 61
  • 2

1 Answers1

3

One option is to use arrayfun.

sig=1;

%sample 3x3 matrix
R=magic(3)/10;

[x,y,z] = ndgrid(-2:2,-2:2,-2:2);
kernel = arrayfun(@(x, y, z) exp(-(norm([x,y,z]*R/sig)^2)/2), x,y,z);

Explanation:

arrayfun takes a function that acts on scalar input(s) to produce scalar output(s), as well as arrays of inputs to pass to the function. Then it loops over input arrays, runs your function on each one, and puts the output of each in the corresponding entry in the output matrix(-es). So arrayfun basically does the nested looping that's slowing everything down for you.

In this example, I also used a anonymous function (also called lambda function) to do the work in the inner-most loop. Since lambda functions need to be a single expression in Matlab, I had to rewrite the inner loop a little bit (using simple algebraic manipulations). If you need to use arrayfun with a function that isn't easily expressed as a lambda, you can always write that function in a separate .m file and pass it to arrayfun.

EDIT: note that you don't have to pre-allocate kernel anymore, and you also don't need offset.

Nicu Stiurca
  • 8,747
  • 8
  • 40
  • 48
  • Great answer, thanks! Works like a charm. What's strange though, the loop works several times faster in this case... Apparently arrayfun is not optimized too well. – user3299285 Feb 12 '14 at 20:21
  • @user3299285 For me, the arrayfun version was much faster. A lot depends on how many loop iterations you have, and how much work you do in each iteration. For such a small loop, (27 iterations total), it's very fast to begin with in absolute terms, so there is not much to gain. Also, norm takes the square root (which is slowish compared to other arithmetic ops) only to come right back and square the results right afterwards. If you want a better comparison of the speed of loop vs arrayfun, replace the body of the loop with the same function used by arrayfun. – Nicu Stiurca Feb 12 '14 at 21:01
  • 1
    I'm afraid I can't agree. Even with replacing the body with the same function I got better results than with arrayfun. It seems that arrayfun is indeed not too fast: http://stackoverflow.com/questions/12522888/arrayfun-can-be-significantly-slower-than-an-explicit-loop-in-matlab-why – user3299285 Feb 13 '14 at 20:13
  • @user3299285 Wow, I guess arrayfun isn't a good solution after all. – Nicu Stiurca Feb 13 '14 at 22:44