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.

Graphics

Graphics in MATLAB is well illustrated in Chapters 12 and 13 of Kevin Sheppard's doc. I leave that to you. If you have any queries or questions even after that, let me know. All further assignments will assume familiarity with basic plotting functions. 

M - Files: III

Debugging

There are three ways of debugging a user-defined function:

1. Brute Force 1: Break the whole function into its (logical) parts and examine what happens as you change inputs and go through a set of commands one by one. May be just copy and paste the sequence of instructions on the command window and what happens to the values of the variables created / modified within the function. It's equivalent to writing the function as a script or breaking it down to smaller scripts.

2. Brute Force 2: Just remove the semicolon at the end of all commands in the function. This will print on the command screen anything or everything that your function does to the variables defined / created within its body. It can potentially, however, result in a very cluttered command window making it difficult to figure out where the program is going wrong (especially when the number of commands is large and you are working with empirical data).

3. Use keyboard & return: Use of the command keyboard within an M-file interrupts the execution of the function and transfers control to the user. In that case, the command prompt symbol is changed to 'K>>'. One can examine or change variables (all MATLAB commands remain valid). To terminate the keyboard mode use the command return.

A related command is pause. You have to enable it first though (use pause on). When an M-file encounters pause it halts execution temporarily until the user presses any key. This can also be a useful debugging tool - or may be just to see after a set of steps if your function is doing what it is supposed to do.

Requesting User Input

It's often needed (or efficient or just plain simple useful) to be able to ask user for inputs (as to whether to continue beyond a point or what values to give to any of the variables or where to read the data from or what name to use to save a variable or the workspace etc etc). The syntax is:

>> user_entry = input('prompt')

A couple of useful extensions include:

>> user_entry = input('prompt', 's')

% In this case the user_entry is taken as a string (even if the user input is numeric)

>> user_entry = input('prompt\n')

% In this case the user_entry input is taken on the 'next line'. Try yourself to see the difference.

An Example

1. Solving the equation f (x) = exp (- x) - sin (x)

function [x, n] = newton1(x0, tol, maxit)

% Using the Newton’s method for solving the equation:

% f (x) = exp (-x) - sin (x)

x = x0; n = 0; done = 0;

while ~done

n = n + 1;

x_new = x -(exp(-x) - sin(x)) / (-exp(-x) - cos(x));

done = (n >= maxit) ( lt (abs(x_new - x), tol);

x = x_new;

end

2. The command feval

The command feval executes the function at pre-specified (argument) values.

For example say you have a function called myfun which takes a single argument 'x'. Then

>> feval('myfun', x0)

evaluates myfun(x) at value x0.

See help feval for more details.

3. Let our myfun be a function that evaluates the function and also calculates the first derivative as in newton1 above:

function [fx, fx_prime] = myfun(x)

% Evaluate f (x) = exp (x) - sin(x) and its first derivative

fx = exp (-x) - sin (x);

fx_prime = - exp (- x) - cos (x);

4. Now let's use myfun(x) and feval together to modify the newton1 function above.

(This also illustrates that you can use a custom function created as an M-file (in this case myfun(x)) inside another custom function (in this case newton2). It's not only more efficient and neater that way, it also makes for easier debugging. And of course, it also illustrates the use of feval)

function [x, n] = newton2 (fname, x0, tol, maxit)

% Using the Newton’s method for solving the equation/function given by fname

% function fname to be entered as a string

% For example if the function name fname, is myfun, it should be entered as 'myfun'

x = x0; n = 0; done = 0;

while ~done

n = n + 1;

[fx, fx_prime] = feval ('myfun', x);

x_new = x – fx / fx_prime;

done = (n > maxit) ( lt (abs(x_new - x), tol);

x = x_new;

end

Note that this now becomes a general function for solving any equation by Newton's method - of course, as long as the function fname is appropriately created. To the extent possible, this is how you should try to create your own functions. Divide the whole into logically sensible parts, and make them as general as possible.