Large scale numerical calculations can put heavy demands on your computer. Some software companies argue that they do not have to do a good job of performance tuning because computer hardware is advancing so rapidly. While it may be possible to foist slug-like word-processors onto the average computer user, engineers and scientists are usually less tolerant of poorly designed and implemented software.
In technical markets the gains in computing performance raise, not lower, the expectations of users. It is often the case that improved hardware and software allows simulations to be run with more detail, which often translates to more realistic results.
MATLAB programs are interpretted. This would seem to make it inapproapriate for large scale scientific computing. The power of MATLAB is realized with its extensive set of libraries which are compiled or are carefully coded in MATLAB to utilize ``vectorization''. The concept of vectorization is central to understanding how to write efficient MATLAB code.
Vectorized code takes advantage, wherever possible, of operations involving data stored as vectors. This even applies to matrices since a MATLAB matrix is stored (by columns) in contiguous locations in the computer's RAM. The speed of a numerical algorithm in MATLAB is very sensitive to whether or not vectorized operations are used.
This section presents some basic considerations to writing efficient MATLAB routines. It is possible, of course, to become obsessed with performance and waste man-hours of programming time in order to save a few seconds of execution time. The ideas presented below require little extra effort, yet can yield significant speed improvements. Once you understand how to use these procedures they will become natural parts of your MATLAB programming style.
For a more advanced discussion of vectorization techniques be sure to check the Mathworks Home Page and in particular their technical notes
The following topics are covered in this section
- Using vector operations instead of loops
- Pre-allocating memory for vectors and matrices
- Correct code is most important
Using vector operations instead of loops
Consider the following loop, translated directly from Fortran or C
dx = pi/30;nx = 1 + 2*pi/dx;for i = 1:nxx(i) = (i-1)*dx;y(i) = sin(3*x(i));end
The preceding statements are perfectly legal MATLAB statements, but they are an inefficient way to create the x and y vectors. Recall that MATLAB allocates memory for variables on the fly. On the first time through the loop (i=1), MATLAB creates two row vectors x and y, each of length one. On each subsequent pass through the loop MATLAB appends new elements to x and y. Not only does this incur extra overhead in the memory allocation calls, the elements of x and y will not be stored in contiguous locations in RAM. Thus, any subsequent operations with x and y, even though these operations may be vectorized, will take a performance hit because of memory access overhead.
The preferred way to create the same two x and y vectors is with the following statements.
x = 0:pi/30:2*piy = sin(3*x);
The first statement creates the row vector, x, with 1 + pi/15 elements stored in contiguous locations in RAM. The second statement creates a new (matrix) variable, y, with the same number of rows and columns as x. Since x is a row vector, as determined by the preceding step, y is also a row vector. If x were, for example, a 5 by 3 matrix, then y = sin(3*x) would create a 5 by 3 matrix, y.
MATLAB is designed to perform vector and matrix operations efficiently. To take maximum advantage of the computer hardware at your disposal, therefore, you should use vectorized operations as much as possible.
Pre-allocating memory for vectors and matrices
Though MATLAB will automatically adjust the size of a matrix (or vector) it is usually a good idea to pre-allocate the matrix. Pre-allocation incurs the cost of memory allocation just once, and it guarantees that matrix elements will be stored in contiguous locations in RAM (by columns).
Consider the following (admitedly artificial) sequence of statements.
dx = pi/30;nx = 1 + 2*pi/dx;nx2 = nx/2;nx = 1 + 2*pi/dx;nx2 = nx/2;for i = 1:nx2x(i) = (i-1)*dx;y(i) = sin(3*x(i));endfor i = nx2+1:nxx(i) = (i-1)*dx;y(i) = sin(5*x(i));end
Since we know the size of x and y, a priori, we can pre-allocate memory for these vectors. Pre-allocation involves creating a matrix (or vector) with one vectorized statement before any of the matrix elements are referenced individually. The ones and zeros functions are typically used to pre-allocate memory.
Here is an improvement of the preceding statements with pre-allocation of memory for x and y.
dx = pi/30;nx = 1 + 2*pi/dx;nx2 = nx/2;for i = 1:nx2nx = 1 + 2*pi/dx;nx2 = nx/2;x = zeros(1,nx); % pre-allocate row-vectors, xy = zeros(1,nx); % and yfor i = 1:nx2x(i) = (i-1)*dx;y(i) = sin(3*x(i));endfor i = nx2+1:nxx(i) = (i-1)*dx;y(i) = sin(5*x(i));end
The statements x(i) = ..., and y(i) = ... still do not take advantage of possibilities for vectorization, but at least the elements of x and y are stored contiguously in RAM. We will improve the efficiency of the loops shortly. First, however, note that we could have written the pre-allocation statements as
x = zeros(1,nx); % pre-allocate row-vectors, xy = x; % and y
The statement y = x does not mean that y will stay equal to x. It simply creates another matrix y with the same ``shape'' as x. Understanding that pre-allocation is important for efficiency will help you understand these apparently confusing twists of MATLAB programming logic.
We can further improve our loops by pulling the assignment of x out of the loops.
x = 0:pi/30:2*pi; % vectorized calculation of xnx = length(x);nx2 = nx/2;for i = 1:nx2nx = length(x);nx2 = nx/2;y = x; % pre-allocate memory for yfor i = 1:nx2y(i) = sin(3*x(i));endfor i = nx2+1:nxy(i) = sin(5*x(i));end
Finally, if we're obsessed with performance, we observe that the calculation of y can also be vectorized. To do this we use the colon notation to refer to segments of the x and y vectors.
x = 0:pi/30:2*pi; % vectorized calculation of xnx = length(x);nx2 = nx/2;nx2 = nx/2;Y (1: nx2) = sin (3*x (1: nx2)); % compute first part of ynx = length(x);nx2 = nx/2;y = x; % pre-allocate memory for yy(1:nx2) = sin(3*x(1:nx2)); % compute first part of yy(nx2+1:nx) = sin(5*x(nx2+1:nx)); % and the second part
To those new to MATLAB programming, the preceding statements may appear unecessarily obfuscated. The comment statements, of course, help, but the logic behind the logic comes from a true understanding of vectorization. Once you get the hang of MATLAB's colon notation you too will come to write code like this. Whenever the speed of MATLAB code is important, there is no substitute for vectorization.
