Sunday, November 16, 2008

MATLAB Tips and Tricks: I


Vectorization

Now would be a good time to talk about how we can optimize performance and write faster codes in MATLAB. The trick is to avoid working on individual elements (as in a for loop) and work on arrays/vectors to the extent possible. This is called vectorization. 

Now while the principle is easily stated, in practice whether a part of the code can be vectorized or not depends on the problem. Even the 'how to' part depends on the nature of the problem. There can be different ways for the same problem and one technique may be better than the other. Vectorization is a skill that develops over time. It helps to be aware that most iterative operations can be (fully or partially) vectorized. 

There are many utilities in MATLAB that makes this task easier (for example, find, end, the colon operator ':', linspace, min, max, cumsum, cumprod etc). The best way to see this is through examples.

Examples

1. Say you have a 3 X 3 matrix and you want to make any value greater than a particular number (say, 4) to be zero. Let:

>> x = magic(3);

% creates a magic square of dimension 3 X 3 which gives:

x =
8 1 6
3 5 7
4 9 2

Now there are two ways of achieving this objective

a) the 'usual' way:

for i = 1 : 3
for j = 1 : 3,
if (x(i, j) > 4),
x(i, j) = - x (i, j);
end
end
end

b) the vectorized version:

index = find(x > 4);
x (index)= - x (index);

The difference is obvious (it's much faster too).

2. Given a set of points, calculate the minimum distance from the origin:

a) without vectorization

function d = minDistance(x, y, z)

% Given a set of of 3-D points specified by columns x, y and z
% compute the minimum distance to the origin

nPoints = length(x);
d = zeros(nPoints, 1);

for k = 1 : nPoints
d(k) = sqrt (x (k)^2 + y (k) ^2 + z (k) ^ 2);
end

d = min(d);

b) with vectorization

function d = minDistance(x, y, z)

% Given a set of of 3-D points specified by columns x, y and z
% compute the minimum distance to the origin

d = sqrt (x .^ 2 + y .^ 2 + z .^ 2);
d = min(d);

3. Tony's Trick: To duplicate a vector of size m X 1 and create a matrix of size m X n

Let 

>> a = [1 2 3];
>> b = a ([1 1 1], :);

gives

>> b = 
1 2 3
1 2 3
1 2 3

Of course, this could have been accomplished by

>> b = a * ones (1, 3)

too. But this requires an unnecessary (and slower) matrix mutliplication.

What is happening here is that the first row of the vector is being indexed n (in this case, 3) times.  In another example, we look at how to manipulate a scalar

>> a = 4;
>> a = a (ones (3, 3));

gives

>> a = 
4 4 4
4 4 4
4 4 4

Again, one avoids the use of matrix multiplication. Similar task, of course, can be achieved by the function repmat too (there is no difference in speed - in fact, repmat employs the Tony's Trick).

I'll stop here on this. You may want to bookmark the Mathworks page on vectorization. Time permitting, we'll cover some more examples of vectorization and, may be, pre-allocation too.

No comments: