Writing Fast MATLAB Code
Pascal Getreuer, June 2006
Contents
1 The Profiler 2
2 Array Preallocation 3
3 Vectorization 5
3.1 Vectorized Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 Vectorized Logic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
4 Inlining Simple Functions 10
5 Numerical Integration 12
5.1 One-Dimensional Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5.2 Multidimensional Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
6 Signal Processing 16
7 Referencing Operations 19
8 Miscellaneous Tricks 22
8.1 Clip a value without using if statements . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
8.2 Convert any array into a column vector . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
8.3 Find the min/max of a matrix or N-d array . . . . . . . . . . . . . . . . . . . . . . . . . 22
8.4 Vector Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
8.5 Flood filling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
8.6 Vectorized use of set on GUI objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
9 Further Reading 24
Introduction
The Matlab programming language is parsed; code is interpreted during runtime. Languages like
C and Fortran are faster because they are first compiled into the computer’s native language. The
advantages of parsing in realtime are greater platform independence, greater language flexibility, and
easier debugging. The disadvantages are slower speed, increased overhead, and limited low-level control.
This article discusses strategies for improving the speed of Matlab code. Keep in mind that Matlab
has gone through many versions and that it is available on many platforms. The fastest method on one
system may not be the fastest on another. This article provides methods that are generally fast, but
makes no claim on what is fastest.
Caution!
Learn the language first. Optimization requires comfort with the syntax and functions of the
language. This article is not a tutorial on MATLAB.
Use comments. Optimized code tends to be terse and cryptic. Help others and yourself by
remembering to comment.
Don’t optimize code before its time. Is optimization worth the effort? If the code will soon
be revised or extended, it will be rewritten anyway.
Only optimize where necessary. Focus your efforts on bottlenecks.
1
1 The Profiler
Matlab 5.0 and newer versions include a tool called the “profiler” that helps determine where the
bottlenecks are in a program. Consider the following function:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
function result = example1(Count)
for k = 1:Count
result(k) = sin(k/50);
if result(k) < −0.9
result(k) = gammaln(k);
end
end
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
To analyze the efficiency this function, first enable the profiler and clear any old profiler data:
>> profile on
>> profile clear
Now run the program. Change the input argument higher or lower so that it takes about a second.
>> example1(5000);
Then enter
>> profile report
The profiler generates an HTML report on the function and launches a browser window. Depending on
the system, profiler results may be a little different from this example.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
MATLAB Profile Report: Summary
Report generated 30-Jul-2004 16:57:01
Total recorded time: 3.09 s
Number of M-functions: 4
Clock precision: 0.016 s
Function List
Name Time Time Time/call Self time Location
example1 3.09 100.0% 1 3.094000 2.36 76.3% ~/example1.m
gammaln 0.73 23.7% 3562 0.000206 0.73 23.7% ../toolbox/matlab/specfun/gammaln.m
profile 0.00 0.0% 1 0.000000 0.00 0.0% ../toolbox/matlab/general/profile.m
profreport 0.00 0.0% 1 0.000000 0.00 0.0% ../toolbox/matlab/general/profreport.m
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
Clicking the “example1” link gives more details:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
100% of the total time in this function was spent on the following lines:
3: for k = 1:Count
2.11 68% 4: result(k) = sin(k/50);
5:
0.14 5% 6: if result(k) < -0.9
0.84 27% 7: result(k) = gammaln(k);
8: end
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The most time-consuming lines are displayed, along with time, time percentage, and line number. The
most costly lines are the computations on lines 4 and 7.
2 Array Preallocation
Matlab’s matrix variables have the ability to dynamically augment rows and columns. For example,
>> a = 2
a =
2
>> a(2,6) = 1
a =
2 0 0 0 0 0
0 0 0 0 0 1
Matlab automatically resizes the matrix. Internally, the matrix data memory must be reallocated with
larger size. If a matrix is resized repeatedly—like within a for loop—this overhead can be significant.
To avoid frequent reallocations, “preallocate” the matrix with the zeros command. Consider the code:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
a(1) = 1;
b(1) = 0;
for k = 2:8000
a(k) = 0.99803 * a(k − 1) − 0.06279 * b(k − 1);
b(k) = 0.06279 * a(k − 1) + 0.99803 * b(k − 1);
end
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The profiler timed this code to take 0.47 seconds. After the for loop, both arrays are row vectors of
length 8000, thus to preallocate, create empty a and b row vectors each with 8000 elements.
3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
a = zeros(1,8000); % Preallocation
b = zeros(1,8000);
a(1) = 1;
b(1) = 0;
for k = 2:8000
a(k) = 0.99803 * a(k − 1) − 0.06279 * b(k − 1);
b(k) = 0.06279 * a(k − 1) + 0.99803 * b(k − 1);
end
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
With this modification, the code takes only 0.14 seconds (over three times faster). Preallocation is often
easy to do, in this case it was only necessary to determine the right preallocation size and add two lines.
What if the final array size can vary? Then use the upper bound on the array size and cut the excess
after the loop:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
a = zeros(1,10000); % Preallocate
count = 0;
for k = 1:10000
v = exp(rand(1)*rand(1));
if v > 0.5 % Conditionally add to array
count = count + 1;
a(count) = v;
end
end
a = a(1:count); % Trim the result
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The average run time of this program is 0.42 seconds without preallocation and 0.18 seconds with it.
Preallocation can also be done with cell arrays, using the cell command to create the desired size.
Using preallocation with a frequently resizing cell array is even more beneficial than with double arrays.
4
3 Vectorization
To “vectorize” a computation means to replace parallel operations with vector operations. This strategy
often improves speed ten-fold. Good vectorization is a skill that must be developed; it requires comfort
with Matlab’s language and creativity.
3.1 Vectorized Computations
Many standardMatlab functions are “vectorized,” they can operate on an array as if the function had
been applied individually to every element.
>> sqrt([1,4;9,16])
ans =
1 2
3 4
>> abs([0,1,2,−5,−6,−7])
ans =
0 1 2 5 6 7
Consider the following function:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
function d = minDistance(x,y,z)
% Find the min distance between a set of points and the origin
nPoints = length(x);
d = zeros(nPoints,1); % Preallocate
for k = 1:nPoints % Compute distance for every point
d(k) = sqrt(x(k)ˆ2 + y(k)ˆ2 + z(k)ˆ2);
end
d = min(d); % Get the minimum distance
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
For every point, the distance between it and the origin is computed and stored in d. The minimum
distance is then found with min. To vectorize the distance computation, replace the for loop with
vector operations:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
function d = minDistance(x,y,z)
% Find the min distance between a set of points and the origin
d = sqrt(x.ˆ2 + y.ˆ2 + z.ˆ2); % Compute distance for every point
d = min(d); % Get the minimum distance
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The modified code performs the distance computation with vector operations. The x, y and z arrays
are first squared using the per-element power operator, .^ (the per-element operators for multiplication
and division are .* and ./). The squared components are added with vector addition. Finally, the
square root of the vector sum is computed per element, yielding an array of distances.
5
The first version of the minDistance program takes 0.73 seconds on 50000 points. The vectorized
version takes less than 0.04 seconds, more than 18 times faster.
Some useful functions for vectorizing computations:
min, max, repmat, meshgrid, sum, cumsum, diff, prod, cumprod, filter
3.2 Vectorized Logic
The previous section shows how to vectorize pure computation. Bottleneck code often involves condi-
tional logic. Like computations, Matlab’s logic operators are vectorized:
>> [1,5,3] < [2,2,4]
ans =
1 0 1
Two arrays are compared per-element. Logic operations return “logical” arrays with binary elements.
How is this useful? Matlab has a few powerful functions for operating on logical arrays:
find: Find indices of nonzero elements.
any: True if any element of a vector is nonzero (or per-column for a matrix).
all: True if all elements of a vector are nonzero (or per-column for a matrix).
>> find([1,5,3] < [2,2,4])
ans =
1 3
>> find(eye(3) == 1)
ans =
1
5
9
The find function returns the indices where the vector logic operation returns true. In the first example,
1 < 2 is true, 5 < 2 is false, and 3 < 4 is true, so find reports that the first and third comparisons are
true. In the second example, find returns the indices where the identity matrix is equal to one. The
indices 1, 5, and 9 correspond to the diagonal of a 3 by 3 matrix.
The any and all functions are simple but occasionally very useful. For example, any(x(:) < 0)
returns true if any element of x is negative.
6
Example 1: Removing elements
The situation often arises where array elements must be removed on some per-element condition. For
example, this code removes all NaN and infinite elements from an array x:
i = find(isnan(x) | isinf(x)); % Find bad elements in x
x(i) = []; % and delete them
Alternatively,
i = find(˜isnan(x) & ˜isinf(x)); % Find elements that are not NaN and not infinite
x = x(i); % Keep those elements
Both of these solutions can be further streamlined by using logical indexing:
x(isnan(x) | isinf(x)) = []; % Delete bad elements
or
x = x(˜isnan(x) & ˜isinf(x)); % Keep good elements
Example 2: Piecewise functions
The sinc function has a piecewise definition,
sinc(x) =
{
sin(x)/x, x 6= 0
1, x = 0
This code uses find with vectorized computation to handle the two cases separately:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
function y = sinc(x)
% Computes the sinc function per−element for a set of x values.
y = ones(size(x)); % Set y to all ones, sinc(0) = 1
i = find(x ˜= 0); % Find nonzero x values
y(i) = sin(x(i)) ./ x(i); % Compute sinc where x ˜= 0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
An interesting alternative is y = (sin(x) + (x == 0))./(x + (x == 0)).
Example 3: Drawing images with meshgrid
The meshgrid function takes two input vectors and converts them to matrices by replicating the first
over the rows and the second over the columns.
>> [x,y] = meshgrid(1:5,1:3)
x =
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
7
y =
1 1 1 1 1
2 2 2 2 2
3 3 3 3 3
The matrices above work like a map for a width 5, height 3 image. For each pixel, the x-location can
be read from x and the y-location from y. This may seem like a gratuitous use of memory as x and y
simply record the column and row positions, but this is useful. For example, to draw an ellipse,
% Create x and y for a width 150, height 100 image
[x,y] = meshgrid(1:150,1:100);
% Ellipse with origin (60,50) of size 15 x 40
Img = sqrt(((x−60).ˆ2 / 15ˆ2) + ((y−50).ˆ2 / 40ˆ2)) > 1;
% Plot the image
imagesc(Img); colormap(copper);
axis image; axis off;
Drawing lines is almost the same, just a change in the formula.
[x,y] = meshgrid(1:150,1:100);
% The line y = x*0.8 + 20
Img = (abs((x*0.8 + 20) − y) > 1);
imagesc(Img); colormap(copper);
axis image; axis off;
Polar functions can be drawn by first converting x and y variables with the cart2pol function.
[x,y] = meshgrid(1:150,1:100);
[th,r] = cart2pol(x − 75,y − 50); % Convert to polar
% Spiral centered at (75,50)
Img = sin(r/3 + th);
imagesc(Img); colormap(hot);
axis image; axis off;
8
Example 4: Polynomial interpolation
Given n points x1, x2, x3, . . . xn and n corresponding function values y1, y2, y3, . . . yn, the coefficients
c0, c1, c2, . . . cn−1 of the interpolating polynomial of degree n− 1 can be found by solving
x1
n−1 x1n−2 · · · x12 x1 1
x2
n−1 x2n−2 · · · x22 x2 1
...
...
xn
n−1 xnn−2 · · · xn2 xn 1
cn−1
cn−2
...
c0
=
y1
y2
...
yn
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
function c = polyint(x,y)
% Given a set of points and function values x and y,
% computes the interpolating polynomial.
x = x(:); % Make sure x and y are both column vectors
y = y(:);
n = length(x); % n = Number of points
%%% Construct the left−hand side matrix %%%
xMatrix = repmat(x, 1, n); % Make an n by n matrix with x on every column
powMatrix = repmat(n−1:−1:0, n, 1); % Make another n by n matrix of exponents
A = xMatrix .ˆ powMatrix; % Compute the powers
c = A\y; % Solve matrix equation for coefficients
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The strategy to construct the left-hand side matrix is to first make two n by n matrices of bases and
exponents and then put them together using the per element power operator, .^ . The repmat function
(“replicate matrix”) is used to make the base matrix xMatrix and the exponent matrix powMatrix.
xMatrix =
x(1) x(1) · · · x(1)
x(2) x(2) · · · x(2)
...
...
x(n) x(n) · · · x(n)
powMatrix =
n− 1 n− 2 · · · 0
n− 1 n− 2 · · · 0
...
...
n− 1 n− 2 · · · 0
The xMatrix is made by repeating the column vector x over the columns n times. The powMatrix is a
row vector with elements n − 1, n − 2, n − 3, . . . , 0 repeated down the rows n times. The two matrices
could also have been created with [powMatrix, xMatrix] = meshgrid(n-1:-1:0, x).
This function is only an example; use the standard polyfit function for serious polynomial interpola-
tion. It is more general and algorithmically more efficient (see polyfit.m).
9
4 Inlining Simple Functions
Every time an M-file function is called, the Matlab interpreter incurs some overhead. Additionally,
many M-file functions begin with conditional code that checks the input arguments for errors or deter-
mines the mode of operation.
Of course, this overhead is negligible for a single function call. It should only be considered when the
function being called is an M-file, the function itself is “simple,” that is, implemented with only a few
lines, and called frequently from within a loop.
For example, this code calls the M-file function median repeatedly to perform median filtering:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% Apply the median filter of size 5 to signal x
y = zeros(size(x)); % Preallocate
for k = 3:length(x)−2
y(k) = median(x(k−2:k+2));
end
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Given a 2500-sample array for x, the overall run time is 0.42 seconds.
“Inlining a function” means to replace a call to the function with the function code itself. Beware that
inlining should not be confused with Matlab’s “inline” function datatype.
Studying median.m (type “edit median” on the console) reveals that most of the work is done using
the built-in sort function. The median call can be inlined as
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% Apply the median filter of size 5 to signal x
y = zeros(size(x)); % Preallocate
for k = 3:length(x)−2
tmp = sort(x(k−2:k+2));
y(k) = tmp(3); % y(k) = median(x(k−2:k+2));
end
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Now the overall run time for a 2500-sample input is 0.047 seconds, nearly 9 ti
本文档为【编写快速的MATLAB代码】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑,
图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。