CODEWITHSUNDEEP
performancematlabvectorvectorization
user3299285
user3299285

Reputation: 61

How to vectorize the 'for' loop in Matlab

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!

Upvotes: 2

Views: 591

Answers (1)

Nicu Stiurca
Nicu Stiurca

Reputation: 8677

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.

Upvotes: 3

Related Questions