Flow of Control and Vectorization in Matlab
by Matt Dunham
Matlab’s constructs for controlling the conditional execution of commands are very similar to those found in most popular programming languages. One major difference, however, is that variables created within the if, for, while, switch, and try statements, are not locally scoped but instead share their scope with all variables in the same function. This is quite unlike java, for example, where a variable created inside a loop can only be used inside that loop.
Contents
- Useful Functions
- if, else, elseif
- for loops
- while loops
- progress and message bars
- switch statements
- try/catch statements
- preallocation
- vectorization
- vectorization examples
- more tips
Useful Functions
zeros, continue, break, return waitbar, inputdlg, listdlg warning, error, rethrow tic, toc bsxfun, repmat, cellfun, cell2mat, mat2cell, arrayfun, vectorize, meshgrid, blkdiag, cumsum, cumprod, filter, conv, accumarray
if, else, elseif
Matlab if statements allow you to execute different code depending on the current state of the program, i.e. the values of certain variables.
test1 = true; test2 = false; test3 = false;
if(test1), A = 1; end % simple if statement on one line. A=1 executed if test1 is true
if(test1) A = 2; % executed if test1 = true else A = 3; % executed if test1 = false end
if(test1) A = 3; % executed if test1 = true elseif(test2) A = 4; % executed if test1 = false and test2 = true elseif(test3) A = 5; % executed if test1 = false, test2 = false, test3 = true else A = 6; % executed if test1=test2=test3=false end
All if statements must end with an end statement.
for loops
For loops allow you to execute a block of code a specified number of times. That number can be determined dynamically as the program runs.
n = ceil(100*rand); % can be set dynamically A = zeros(n,1); % improve speed by preallocating space for i=1:n % set i = 1, then loop and increment i by 1, until i = n A(i,1) = max(i,50); % execute code within the loop - usually depends on i. end % both i and A can then be accessed outside the loop.
We can nest for loops and include if statements. Note, we do not have to start looping from 1. We feature the continue, break, and return commands. Continue instructs Matlab to skip directly to the beginning of the current loop without first executing the lines directly below it. Break, on the other hand, breaks from the current loop completely. Return breaks completely from the current script or function without executing any further code.
A = rand(20,20,20); counter = 0; for i=1:size(A,1) for j=2:size(A,2) for k=3:size(A,3) %if k is even, go immediately to beginning of loop if(mod(k,2) == 0), continue; end %if j+k is prime, break from this inner loop completely if(isprime(j+k)), break ; end %if all three of i,j,k prime, stop all further execution. if(all(isprime([i,j,k])) && false), return ; end if(isprime(floor(100*A(i,j,k)))) counter = counter + 1; end end end end
The continue, break, and return statements should be used sparingly as they can easily obscure the code and can almost always be replaced by if,else,elseif statements.
while loops
While loops are used to execute a block of code until some condition is satisfied. This condition is usually more complicated than simply reaching a set number of iterations as with a for loop. The comments on scope, and the continue, break and return statements apply equally to while loops.
A = true; B = true; C = true; val = 1; while(A || B || C) val = 2*val +1; A = isprime(val); B = val < 10; C = ((round(sqrt(val)))^2) == val; end
Here is common code idiom involving break. This effectively allows us to test at the end of the loop, (or in multiple spots).
while(true) %execute code if(condition) break; end end
progress and message bars
Newer versions of Matlab have nice, easy to use, graphical message and progress windows. Search for Predefined Dialog Boxes in help for a full list of available dialogs. We give examples of three here, waitbar(), inputdlg()_, and listdlg().
Keep the user informed as to the progress of a loop or lengthy calculation.
w = waitbar(0,'My Progress Bar'); % create a new waitbar, w with 0% progress for i=1:500 isprime(i); w = waitbar(i/500,w,['iteration: ',num2str(i)]); % update the wait bar each iteration end close(w); % remember to close it
Ask the user for input with a graphical window
answer = inputdlg('What would you like for dinner?'); % ask the user for input
Give the user a list of options. The index of the selected one is returned.
options = {'Chicken','Beef','Fish','Pasta'}; message = 'Here are your options:'; [selection, ok] = listdlg('PromptString',message,'ListString',options);
switch statements
Switch statements are useful when what code to execute depends on a variable that takes on a countable number of values. Most commonly, this value is an integer or a string. Switch statements can be replaced by a long series of if-else statements but this usually results in less readable code. Note that unlike languages such as C or java, switch statements do not fall through; that is, the code from, (at most), one case statement is executed. As such, break statements are not necessary.
color = 'blue'; switch color % switching variable case 'red' A = 1; % code for case 'red' case 'blue' A = 2; case {'green','purple'} % either 'green' or 'purple' A = 3; otherwise % optional 'catch all' A = 4; end
try/catch statements
Try catch blocks give you some control over Matlab error handling. They are useful for executing code that might potentially fail, such as writing to a file, allowing you to perform cleanup or recover gracefully. Overuse, however, can easily obscure your program and does not necessarily make it easier to debug.
a = rand; b = a*(a< 0.5); try c = a / b; assert(true); % set to false to have code throw an error catch ME % disaster recovery, cleanup, inform user, etc... display('Something went wrong'); warning('WARNING:ID','my own warning message'); display(ME.message); % ME is a structure with info on the error %error('my own error message'); % stops execution %rethrow(ME); % rethrows the original error and stops execution end
preallocation
Matlab stores matrices in contiguous blocks of memory. When the size of a matrix changes, Matlab, if it has not preallocated enough space, must find a new chunk of memory large enough and copy the matrix over. When a matrix grows inside of a loop, this process may have to be repeated over and over again causing huge delays. It can therefore significantly speed up your code by preallocating a chunk of memory before entering into a loop. The zeros() command is the most common way to do this. Below we see two simple loops in which we store the numbers 1 to 30 000. We preallocate only in the second. Timing the two loops with the tic() and toc() commands we see that preallocating in this case speeds up the code by about 30 times. The larger the matrices, the more important this becomes.
tic for i = 1:30000 A(i) = i; end without = toc
without =
1.3971tic B = zeros(30000,1); % Preallocate B with the zeros command. for i = 1:30000 B(i) = i; end with = toc ratio = without / with
with =
0.0471
ratio =
29.6756vectorization
Matlab is an interpreted language, which means that each line of code must be reduced to machine instructions as the program runs, whereas with compiled code, this is done before execution. Moreover, compiled code can often be optimized automatically by the compiler in ways interpreted code cannot. As such, interpreted code is generally slower than an equivalent compiled version.
Many of Matlab’s built in functions, however, particularly those involving matrix operations, are highly optimized and compiled in a low level language like C or Fortran resulting in very fast code. But where Matlab is at a disadvantage, is with regard to loops. Although recent versions have seen a considerable increase in speed, loops are still a major bottleneck.
Thankfully, we can frequently replace loops with matrix operations or calls to fast, built in functions – a process called vectorization. Learning how to do this well is an extremely important skill.
We will see a number of examples comparing vectorized with non-vectorized code. Many of these might seem obvious to seasoned Matlab programmers but we show to them anyway to emphasize two points: the speed increase, and the relative terseness of the vectorized versions. Unfortunately publishing the document skews the timing results but try the code out for yourself. For some of the examples, we see a 30 fold increase in speed.
vectorization examples
A = rand(200,200); % We will use this as our data
Most functions in Matlab are already vectorized, so that to take the log of every number in an array A, for instance, we simply execute B = log(A).
non-vectorized version
tic % time the code Bnv = zeros(size(A)); % We preallocate to level the playing field for i=1:size(A,1) for j=1:size(A,2); Bnv(i,j) = log(A(i,j)); end end nonvec = toc;
vectorized version
tic Bv = log(A); vec = toc;
assert(isequal(Bnv,Bv)); ratio = nonvec / vec;
Matlab supports parallel indexing and assignment so that we can retrieve and assign multiple values at once.
A1 = A; A2 = A;
rsndx = 1:100; csndx = 80:130; rtndx = 101:200; ctndx = 150:200;
non-vectorized version
tic for i=1:numel(rsndx) for j=1:numel(csndx) A1(rsndx(i),csndx(j)) = A1(rtndx(i),ctndx(j)); end end nonvec = toc;
vectorized version
tic A2(rsndx,csndx) = A2(rtndx,ctndx); vec = toc;
ratio = nonvec / vec; assert(isequal(A1,A2));
Here we see the benefit of logical indexing.
non-vectorized version
tic B1 = []; % note, it is difficult to preallocate here counter = 1; for j=1:size(A,2) for i=1:size(A,1) if(A(i,j) < 0.2) B1(counter,1) = A(i,j); counter = counter + 1; end end end nonvec = toc;
vectorized version
tic B2 = A(A < 0.2); vec = toc;
ratio = nonvec / vec; assert(isequal(B1,B2));
Here we perform three tricks at once as it were. Recall that operators such as ^, \, have element-wise equivalents, (e.g. .^), which we can apply to the corresponding elements of two same-sized matrices. Secondly, Matlab performs automatic scalar expansion in expressions like A+1, and thirdly, we can easily multiply two matrices together without loops. Most loops involving patterned additions and multiplications of vector elements can be translated, with a little thought, into equivalent vectorized statements.
non-vectorized version
tic B1 = zeros(size(A)); for i=1:size(A,1) for j=1:size(A,2) T = 0; for k=1:size(A,1) T = T + A(i,k)*A(j,k); end B1(i,j) = T * (A(i,j)/2) + 1; end end nonvec = toc;
vectorized version
tic B2 = ((A*A') .* (A/2)) + 1; vec = toc;
test = mean(abs(B1(:) - B2(:))); % very small differences between B1, & B2 because of numerical error ratio = nonvec / vec;
Recall from the chapter on matrices that we can use repmat() or bsxfun() to perform element-wise operations on non-scalar matrices of different sizes as long as a singleton dimension can be extended to make them the same size. Here we subtract off the mean of the third dimension and leave our ‘non-vectorized’ version at least somewhat vectorized to emphasize the role of bsxfun().
A3d = rand(100,100,100); A1 = A3d; A2 = A3d; A3 = A3d;
non-vectorized version
tic m = mean(A1,3); for i=1:size(A1,3) A1(:,:,i) = A1(:,:,i) - m; end nonvec = toc;
vectorized version
tic A2 = bsxfun(@minus,A2,mean(A2,3)); vec = toc;
We could have also used repmat() as follows, but this requires more memory and is slightly slower.
tic A3 = A3 - repmat(mean(A3,3),[1,1,size(A3,3)]); rep = toc;
assert(isequal(A1,A2,A3));
Here is a real example involving the last two techniques. We are calculating a scalar value lambda given X, y, & W, with the following dimensions.
- X is n-by-d
- y is n-by-1
- W is d-by-1
With the following formula, we can calculate a single lambda value per W vector.
lambda = 2*max(abs(X'*(y-X*W)),[],1);
However, we would like to calculate multiple lambda values given multiple W vectors. We could use a loop, but then this would not be much of an example. Instead, we stack k of the d-by-1 W vectors, column-wise into a d-by-k matrix before hand. XW is then an n-by-k matrix. We use bsxfun() to effectively expand y, which was an n-by-1 vector, into an n-by-k vector with k identical columns and perform the subtraction. Finally we multiply X’ by the resuling n-by-k matrix yielding a d-by-k matrix, and take the maximum along the first dimension yielding the k lambda values we were after. Note that this new version still works even if k=1, that is, even if we are after only one lambda value.
lambdaVals = 2*max(abs(X'*(bsxfun(@minus,y,X*W))),[],1);
Vector equations tend to naturally generalize into matrix equations.
Here is one final example.
Suppose we have a large numeric matrix and we want to apply a function to arbitrary sized blocks of it. That is, we want to partition a matrix of size m-by-n into many smaller matrices of differing sizes, apply a function to each block, and group the results back together. We could extract each block first with a long series of indexing operations and then loop over them all applying the function, but there is better way involving the mat2cell() and cellfun() functions.
We have not discussed cells yet and you can skip this example until after you read about them in the next chapter if you like. However, for now you could just use this example as a template.
A = rand(100,40); % here is our data
Partition the matrix into 12 blocks of different sizes. These blocks are stored in a 4×3 cell array. Notice the sizes of each of the 12 blocks and how we achieved these sizes with the inputs to mat2cell().
groups = mat2cell(A,[10,30,20,40],[5,27,8])
groups =
[10x5 double] [10x27 double] [10x8 double]
[30x5 double] [30x27 double] [30x8 double]
[20x5 double] [20x27 double] [20x8 double]
[40x5 double] [40x27 double] [40x8 double]Create a function to apply to each block; we will choose something simple like replacing each element in a block with the block’s largest value.
f = @(x)repmat(max(x(:)),size(x));
Use the cellfun() function to apply this function to every one of the 12 elements in groups, (i.e. to every matrix block). We have to set ‘UniformOutput’ to false because the sizes of the elements returned by cellfun() will be different.
groupSums = cellfun(f,groups,'UniformOutput',false)
groupSums =
[10x5 double] [10x27 double] [10x8 double]
[30x5 double] [30x27 double] [30x8 double]
[20x5 double] [20x27 double] [20x8 double]
[40x5 double] [40x27 double] [40x8 double]We then convert back to a numeric matrix with the same size as our original matrix A.
B = cell2mat(groupSums);
more tips
In the last example, we used cellfun() function but there is a similar function arrayfun() that applies a function to every element of an array. When other vectorization techniques fail, this can be a better alternative than looping over every element yourself.
Some functions, like mvnpdf() for example, interpret an n-by-d matrix, not as n-times-d elements but as n, d-dimensional vectors. If this is not what we are after, we can convert the matrix into a vector using the (:) operator, pass it to the function, and reshape the output back into the original size with the reshape() function.
The vectorize() function takes in a string or function handle and converts all operators, (e.g ^) to their element-wise equivalents, (e.g. .^). This can be useful when using someone else’s function that was not vectorized to begin with.
Recall from the matrix chapter that there are many functions that will create matrices for such as meshgrid(), or blkdiag(), yet again helping us avoid loops.
When a value vec(i) depends on on entries v(1)…v(i-1) for instance, we can use functions like cumsum(), cumprod(), filter(), conv(), and accumarray(). See their help entries for more information.
Finally, if a loop is essential and proves to be a significant bottleneck in your program, consider compiling it via emlmex, (if possible). Alternatively, write the loop in another language like C or Java and call that code from within your Matlab program. For details, see the chapter on calling external code.
clear all;
Popularity: 1% [?]
Related Posts
Comments
Tell me what you're thinking...
and oh, if you want a pic to show with your comment, go get a gravatar!
Include MATLAB code in your comment by doing the following:
<pre lang="MATLAB">
%insert code here
</pre>
Translator
Follow Us
Video Posts
