pobieranie: Writing fast Matlab code (ang)
Pobierz dokument pdfPrzeglądaj wersję html pliku:
Pascal Getreuer, February 2009
Contents
Page
1 Introduction 2 The Profiler 3 Array Preallocation 4 JIT Acceleration 5 Vectorization 6 Inlining Simple Functions 7 Referencing Operations 8 Solving Ax = b 9 Numerical Integration 10 Signal Processing 11 Miscellaneous 12 Further Reading
1 3 5 7 9 14 16 19 22 26 29 31
1
Introduction
Matlab is a prototyping environment, meaning it focuses on the ease of development with language flexibility, interactive debugging, and other conveniences lacking in performance-oriented languages like C and Fortran. While Matlab may not be as fast as C, there are ways to bring it closer.
Disclaimer This is not a beginner’s tutorial to Matlab, but a tutorial on
performance. The methods discussed here are generally fast, but no claim is made on what is fastest. Use at your own risk.
Before beginning our main topic, let’s step back for a moment and ask ourselves what we really want. Do we want the fastest possible code?—if so, we should switch to a performance language like C. Probably more accurately, we want to spend less time total from developing, debugging, running, and until finally obtaining results. This section gives some tips on working smart, not hard.
M-Lint Messages
In more recent versions of Matlab, the integrated editor automatically gives feedback on possible code errors and suggested optimizations. These messages are helpful in improving a program and in learning more about Matlab.
for k = 1:NumTrials r = rand; x(k) = runsim(r); x might be growing inside a loop. Consider preallocating for speed. end hist(x); % Plot histogram of simulation outputs
Hold the mouse cursor over underlined code to read a message. Or, see all M-Lint messages at once ::::::::::::::: with Tools → M-Lint → Show M-Lint Report.
Organization
Use a separate folder for each project.
A separate folder for each project keeps related things together and simplifies making copies and backups of project files.
Write header comments, especially H1.
The first line of the header comment is called the H1 comment. Type help(cd) to get a listing of all project files and their H1 comments.
1
Save frequent console commands as a script.
If you find yourself repeating certain commands on the console, save them as a script. Less typing means fewer opportunities for typos.
Avoid losing data
Don’t use clear all in a script.
This is an unfortunate common practice—any important variables in the base workspace will be irretrievably lost.
Beware of clobber.
“File clobber” refers to the kind of data loss when a file is accidentally overwritten with another one having the same filename. This phenomenon can occur with variables as well:
>> result = big operation(input1); ... >> result = big operation(input2); % Store the result of a time−consuming operation % Run again with another input
Variable result was clobbered and the first output was lost.
Beware of what can crash MATLAB.
While Matlab is generally reliable, crashes are possible when using third-party MEX functions or extremely memory-intensive operations, for example, with video and very large arrays. Now with good working habits covered, we begin our discussion of writing fast Matlab code. The rest of this article is organized by topic, first on techniques that are useful in general application, and next on specific computational topics (table of contents is on the first page).
2
2
The Profiler
Matlab 5.0 (R10) and newer versions include a tool called the profiler that helps identify bottlenecks in a program. Use it with the profile command: profile profile profile profile on off clear report Turn the profiler on Turn it back off Clear profile statistics View the results from the profiler
For example, consider profiling 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 and clear the profiler, run the function, and then view the profile report:
>> profile on, profile clear >> example1(5000); >> profile report
There is a slight parsing overhead when running code for the first time; run the test code twice and time the second run. The profiler report command shows a report. Depending on the system, profiler results may differ from this example.
MATLAB Profile Report: Summary
Report generated 30-Jul-2004 16:57:01 Total recorded time: Number of M-functions: Clock precision: 3.09 s 4 0.016 s
3
Function List
Name
example1 gammaln profile profreport
Time
3.09 0.73 0.00 0.00 100.0% 23.7% 0.0% 0.0%
Time
1 3562 1 1
Time/call
3.094000 0.000206 0.000000 0.000000
Self time
2.36 0.73 0.00 0.00 76.3% 23.7% 0.0% 0.0%
Location
~/example1.m ../toolbox/matlab/specfun/gammaln.m ../toolbox/matlab/general/profile.m ../toolbox/matlab/general/profreport.m
Clicking the “example1” link gives more details: Lines where the most time was spent
Line Number Code result(k) = sin(k/50); result(k) = gammaln(k); if result(k) < -0.9 Calls Total Time % Time
4 7 6
Totals
5000 721 5000
2.11 0.84 0.14 3.09
s s s s
68% 27% 5% 100%
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. Another helpful section of the profile report is “M-Lint Results,” which gives feedback from the M-Lint code analyzer. Possible errors and suggestions are listed here. M-Lint results
Line number Message ‘result’ might be growing inside a loop. Consider preallocating for speed. ‘result’ might be growing inside a loop. Consider preallocating for speed.
4 7
(Preallocation is discussed in the next section.) The profiler has limited time resolution, so to profile a piece of code that runs too quickly, run the test code multiple times with a loop. Adjust the number of loop iterations so that the time it takes to run is noticeable. More iterations yields better time resolution. The profiler is an essential tool for identifying bottlenecks and per-statement analysis, however, for more accurate timing of a piece of code, use the tic/toc stopwatch timer.
>> tic; example1(5000); toc; Elapsed time is 3.082055 seconds
For serious benchmarking, also close your web browser, anti-virus, and other background processes that may be taking CPU cycles.
4
3
>> a = 2 a = 2
Array Preallocation
Matlab’s matrix variables have the ability to dynamically augment rows and columns. For example,
>> 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 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
This code takes 0.47 seconds to run. 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.
a = zeros(1,8000); b = zeros(1,8000); a(1) = 1; b(1) = 0; % Preallocation
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.
5
What if the final array size can vary? Here is an example:
a = zeros(1,10000); count = 0; for k = 1:10000 v = exp(rand*rand); if v > 0.5 % Conditionally add to array count = count + 1; a(count) = v; end end a = a(1:count); % Trim extra zeros from the results % Preallocate
The average run time of this program is 0.42 seconds without preallocation and 0.18 seconds with it. Preallocation is also beneficial for cell arrays, using the cell command to create a cell array of the desired size.
6
4
JIT Acceleration
Matlab 6.5 (R13) and later feature the Just-In-Time (JIT) Accelerator for improving the speed of M-functions, particularly with loops. By knowing a few things about the accelerator, you can improve its performance. The JIT Accelerator is enabled by default. To disable it, type “feature accel off” in the console, and “feature accel on” to enable it again. As of Matlab R2008b, only a subset of the Matlab language is supported for acceleration. Upon encountering an unsupported feature, acceleration processing falls back to non-accelerated evaluation. Acceleration is most effective when significant contiguous portions of code are supported.
Data types: Code must use supported data types for acceleration: double (both real and complex), logical, char, int8–32, uint8–32. Some struct, cell, classdef, and function handle usage is supported. Sparse arrays are not accelerated. Array shapes: Array shapes of any size with 3 or fewer dimensions are supported. Changing the shape or data type of an array interrupts acceleration. A few limited situations with 4D arrays are accelerated. Function calls: Calls to built-in functions and M-functions are accelerated. Calling MEX functions and Java interrupts acceleration. (See also page 14 on inlining simple functions.) Conditionals and loops: The conditional statements if, elseif, and simple switch statements are supported if the conditional expression evaluates to a scalar. Loops of the form for k=a:b, for k=a:b:c, and while loops are accelerated if all code within the loop is supported.
In-place computation
Introduced in Matlab 7.3 (R2006b), the element-wise operators (+, .*, etc.) and some other functions can be computed in-place. That is, a computation like
x = 5*sqrt(x.ˆ2 + 1);
is handled internally without needing temporary storage for accumulating the result. An M-function can also be computed in-place if its output argument matches one of the input arguments.
x = myfun(x); function x = myfun(x) x = 5*sqrt(x.ˆ2 + 1); return;
To enable in-place computation, the in-place operation must be within an M-function (and for an inplace function, the function itself must be called within an M-function). Currently, there is no support for in-place computation with MEX-functions.
7
Multithreaded Computation
Matlab 7.4 (R2007a) introduced multithreaded computation for multicore and multiprocessor computers. Multithreaded computation accelerates some per-element functions when applied to large arrays (for example, .^, sin, exp) and certain linear algebra functions in the BLAS library. To enable it, select File → Preferences → General → Multithreading and select “Enable multithreaded computation.” Further control over parallel computation is possible with the Parallel Computing Toolbox using parfor and spmd.
JIT-Accelerated Example
For example, the following loop-heavy code is supported for acceleration:
function B = bilateral(A,sd,sr,R) % The bilateral image denoising filter B = zeros(size(A)); for i = 1:size(A,1) for j = 1:size(A,2) zsum = 0; for m = −R:R if i+m >= 1 && i+m <= size(A,1) for n = −R:R if j+n >= 1 && j+n <= size(A,2) z = exp(−(A(i+m,j+n) − A(i,j))ˆ2/(2*sdˆ2))... * exp(−(mˆ2 + nˆ2)/(2*srˆ2)); zsum = zsum + z; B(i,j) = B(i,j) + z*A(i+m,j+n); end end end end B(i,j) = B(i,j)/zsum; end end
For a 128 × 128 input image and R = 3, the run time is 53.3 seconds without acceleration and 0.68 seconds with acceleration.
8
5
Vectorization
A computation is vectorized by taking advantage of vector operations. A variety of programming situations can be vectorized, and often improving speed to 10 times faster or even better. Vectorization is one of the most general and effective techniques for writing fast M-code.
5.1 Vectorized Computations
Many standard Matlab 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 3 >> abs([0,1,2,−5,−6,−7]) ans = 2 4 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, its distance from the origin is computed and stored in d. For speed, array d is preallocated (see Section 3). 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); d = min(d); % Compute distance for every point % 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. (A further improvement: it is equivalent to compute d = sqrt(min(x.^2 + y.^2 + z.^2))).
9
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, accumarray
5.2 Vectorized Logic
The previous section shows how to vectorize pure computation. But bottleneck code often involves conditional 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 values. How is this useful? Matlab has a few powerful functions for operating on logical arrays:
>> find([1,5,3] < [2,2,4]) % Find indices of nonzero elements ans = 1
3 % True if any element of a vector is nonzero % (or per−column for a matrix)
>> any([1,5,3] < [2,2,4]) ans = 1 >> all([1,5,3] < [2,2,4]) ans = 0
% True if all elements of a vector are nonzero % (or per−column for a matrix)
Vectorized logic also works on arrays.
>> find(eye(3) == 1) ans = 1 5 9 % (element locations given as indices) % Find which elements == 1 in the 3x3 identity matrix
By default, find returns element locations as indices (see page 16 for indices vs. subscripts).
10
Example 1: Vector Normalization
To normalize a single vector v to unit length, one can use v = v/norm(v). However, to use norm for set of vectors v(:,1), v(:,2), . . . requires computing v(:,k)/norm(v(:,k)) in a loop. Alternatively,
vMag = sqrt(sum(v.ˆ2)); v = v./vMag(ones(1,size(v,1)),:);
Example 2: 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)); x(i) = []; % Find bad elements in x % and delete them
Alternatively,
i = find(˜isnan(x) & ˜isinf(x)); x = x(i); % Find elements that are not NaN and not infinite % 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 3: Piecewise functions
The sinc function has a piecewise definition, sinc(x) = sin(x)/x, 1, x=0 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)); i = find(x ˜= 0); y(i) = sin(x(i)) ./ x(i); % Set y to all ones, sinc(0) = 1 % Find nonzero x values % Compute sinc where x ˜= 0
A concise alternative is y = (sin(x) + (x == 0))./(x + (x == 0)).
11
Example 4: 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 1 1 2 2 2 3 3 3 4 4 4 5 5 5 y = 1 2 3 1 2 3 1 2 3 1 2 3 1 2 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 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); % Spiral centered at (75,50) Img = sin(r/3 + th); imagesc(Img); colormap(hot); axis image, axis off
% Convert to polar
12
Example 5: Polynomial interpolation
Given n points x1 , . . . xn and corresponding function values y1 , . . . yn , the coefficients c0 , . . . , cn−1 of the interpolating polynomial P (x) = cn−1 xn−1 + · · · + c1 x + c0 can be found by solving x1 n−1 x2 n−1 . . . n−1 xn
x1 n−2 x2 n−2 xn n−2
··· ··· ···
x1 2 x2 2 xn 2
x1 x2 xn
1 1 . . . 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(:); y = y(:); n = length(x); % Make sure x and y are both column vectors % n = Number of points
% Construct the matrix on the left−hand−side 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×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(2) x(2) · · · . . . x(n) x(n) · · · x(1) x(2) . . . x(n) powMatrix = n−1 n−2 n−1 n−2 . . . n−1 n−2 ··· ··· ··· 0 0 . . . 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). The code above is only an example—use the standard polyfit function for serious use.
13
6
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 determines 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 is 0.047 seconds, nearly 9 times faster. Furthermore, by inlining median, it can be specifically tailored to evaluating 5-sample medians. But this is only an example; if the Signal Processing Toolbox is available, y = medfilt1(x,5) is much faster. A surprising number of Matlab’s functions are implemented as M-files, of which many can be inlined in just a few lines. In Matlab 5.3 (R11), the following functions are worth inlining: conv, cross, fft2, fliplr, flipud, ifft, ifft2, ifftn, ind2sub, ismember, linspace, logspace, mean, median, meshgrid, poly, polyval, repmat, roots, rot90, setdiff, setxor, sortrows, std, sub2ind, union, unique, var The list above is for Matlab R11; some of these functions have since become built-in. Use the which command or try “edit function name” to determine whether a function is implemented as an M-file.
14
For example, in Matlab R11 and earlier, ifft was implemented by simply calling fft and conj. If x is a one-dimensional array, y = ifft(x) can be inlined with y = conj(fft(conj(x)))/length(x). Another example: b = unique(a) can be inlined with
b = sort(a(:)); b( b((1:end−1)') == b((2:end)') ) = [];
While repmat has the generality to operate on matrices, it is often only necessary to tile a vector or just a scalar. To repeat a column vector y over the columns n times,
A = y(:,ones(1,n)); % Equivalent to A = repmat(y,1,n);
To repeat a row vector x over the rows m times,
A = x(ones(1,m),:); % Equivalent to A = repmat(x,m,1);
To repeat a scalar s into an m×n matrix,
A = s(ones(m,n)); % Equivalent to A = repmat(s,m,n);
This method avoids the overhead of calling an M-file function. It is never slower than repmat (critics should note that repmat.m itself uses this method to construct mind and nind). For constructing matrices with constant value, there are other efficient methods, for example, s+zeros(m,n). Don’t go overboard. Inlining functions is only beneficial when the function is simple and when it is called frequently. Doing it unnecessarily obfuscates the code.
15
7
Referencing Operations
Referencing in Matlab is varied and powerful enough to deserve a section of discussion. Good understanding of referencing enables vectorizing a broader range of programming situations.
7.1 Subscripts vs. Indices
Subscripts are the most common method used to refer to matrix elements, for example, A(3,9) refers to row 3, column 9. Indices are an alternative referencing method. Consider a 10 × 10 matrix A. Internally, Matlab stores the matrix data linearly as a one-dimensional, 100-element array of data in column-major order. 1 2 3 . . . 10 11 21 12 22 13 23 . . . . . . 20 30 ··· 81 82 83 . . . 90 91 92 93 . . . 100
···
An index refers to an element’s position in this one-dimensional array. Using an index, A(83) also refers to the element on row 3, column 9. It is faster to scan down columns than over rows. Column-major order means that elements along a column are sequential in memory while elements along a row are further apart. Scanning down columns promotes cache efficiency. Conversion between subscripts and indices can be done with the sub2ind and ind2sub functions. However, because these are M-file functions rather than fast built-in operations, it is more efficient to compute conversions directly (also see Section 6). For a two-dimensional matrix A of size M×N, the conversions between subscript (i,j) and index (index) are A(i,j) ↔ A(i + (j-1)*M) A(index) ↔ A(rem(index-1,M)+1, floor((index-1)/M)+1) Indexing extends to N-D matrices as well, with indices increasing first through the columns, then through the rows, through the third dimension, and so on. Subscript notation extends intuitively, A(. . . , dim4, dim3, row, col).
7.2 Vectorized Subscripts
It is useful to work with submatrices rather than individual elements. This is done with a vector of indices or subscripts. If A is a two-dimensional matrix, a vector subscript reference has the syntax A(rowv, colv)
16
where rowv and colv are vectors. Both may be of any length and their elements may be in any order. If either is a matrix, it is reshaped to a vector. There is no difference between using row vectors or column vectors in vector subscripts. Let M = length(rowv) and N = length(colv). side∗ returns a submatrix of size M×N: A(rowv(1), colv(1)) A(rowv(1), colv(2)) A(rowv(2), colv(1)) A(rowv(2), colv(2)) . . . A(rowv(M), colv(1)) A(rowv(M), colv(2)) Then a vector subscripted matrix on the right-hand
A(rowv(1), colv(3)) · · · A(rowv(2), colv(3)) · · · A(rowv(M), colv(3)) · · ·
A(rowv(1), colv(N)) A(rowv(2), colv(N)) . . . A(rowv(M), colv(N))
A vector subscripted matrix can also appear on the left-hand side as A(rowv, colv) = scalar or M×N matrix If any elements in the destination reference are repeated, the source element with the greater index dominates. For example,
>> A = zeros(2); A([1,2],[2,2]) = [1,2;3,4] A = 0 0 2 4
Vector subscript references extend intuitively in higher dimensions.
7.3 Vector Indices
Multiple elements can also be referenced with vector indices. A(indexv) where indexv is an array of indices. On the right-hand side, a vector index reference returns a matrix the same size and shape as indexv. If indexv is a 3 × 4 matrix, then A(indexv) is the 3 × 4 matrix A(indexv(1, 1)) A(indexv(1, 2)) A(indexv(1, 3)) A(indexv(1, 4)) A(indexv) = A(indexv(2, 1)) A(indexv(2, 2)) A(indexv(2, 3)) A(indexv(2, 4)) A(indexv(3, 1)) A(indexv(3, 2)) A(indexv(3, 3)) A(indexv(3, 4)) While vector subscripts are limited to referring to block-shaped submatrices, vector indices can refer to any subset of elements. If a vector index reference is on the left-hand side, the right-hand side must return a matrix of the same size as indexv or a scalar. As with vector subscripts, ambiguous duplicate assignments use the value with greater source index.
∗ The
left- and right-hand sides (LHS and RHS) refer to the two sides of an assignment, “LHS = RHS.”
17
7.4 Reference Wildcards
Using the wildcard, :, in a subscript refers to an entire row or column. For example, A(:,1) refers to every row in column one–the entire first column. This can be combined with vector subscripts, A([2,4],:) refers to the second and fourth rows. When the wildcard is used in a vector index, the entire matrix is referenced. On the right-hand side, this always returns a column vector. A(:) = column vector This is frequently useful: for example, if a function input must be a row vector, the user’s input can be quickly reshaped into row vector with A(:).’ (make a column vector and transpose to a row vector). A(:) on the left-hand side assigns all the elements of A, but does not change its size. For example, A(:) = 8 changes all elements of matrix A to 8.
7.5 Logical Indexing
Given a logical array mask with the same size as A,
A(mask)
refers to all elements of A where the corresponding element of mask is 1 (true). It is equivalent to A(find(mask)). A common usage of logical indexing is to refer to elements based on a per-element condition, for example,
A(abs(A) < 1e−3) = 0
sets all elements with magnitude less than 10−3 to zero. Logical indexing is also useful to select nonrectangular sets of elements, for example,
A(logical(eye(size(A))))
references the diagonal elements of A. Note that for a right-hand side reference, it is faster to use diag(A) to get the diagonal of A.
7.6 Deleting Submatrices with [ ]
Elements in a matrix can be deleted by assigning the empty matrix. For example, A([3,5]) = [ ] deletes the third and fifth element from A. If this is done with index references, the matrix is reshaped into a row vector. It is also possible to delete with subscripts if all but one subscript are the wildcard. A(2,:) deletes the second row. Deletions like A(2,1) = [ ] or even A(2,1:end) = [ ] are illegal. = []
18
8
Solving Ax = b
Many problems involve solving a linear system, that is, solving the Ax = b problem b1 x1 a1,1 a1,2 · · · a1,n a2,1 a2,2 · · · a2,n x2 b2 . . . . = . . .. . . . . . . . . . . . bn xn an,1 an,2 · · · an,n The condition number cond(A) describes how sensitive the problem is to small perturbations. Matlab can efficiently approximate cond(A) with condest. Smaller condition number is better, and in this case A is said to be well-conditioned. An infinite condition number implies that A is singular—the problem cannot be solved, though it may be solved approximately with the pseudoinverse (x = pinv(A)*b). Small- to moderately-sized problems are solved efficiently by the backslash operator:
x = A\b; % Solves A*x = b
We shall focus on large problems.
8.1 Iterative Methods
For large problems, Matlab is well-equipped with iterative solvers: bicg, bicgstab, cgs, gmres, lsqr, minres, pcg, symmlq, qmr. The general calling syntax (except for gmres) is [x,flag,relres,iter,resvec] = method (A,b,tol,maxit,M1,M2,x0) where the required arguments are the matrix A, the right-hand side b, and the solution x. The minres, pcg, and symmlq methods require A to be symmetric (A = A’) and pcg additionally requires positive definiteness (all eigenvalues are positive). The optional input variables are tol maxit M1, M2 x0 flag relres iter resvec accuracy tolerance (default 10−6 ) maximum number of iterations (default 20) the method uses preconditioner M = M1*M2; it solves M −1 Ax = M −1 b initial guess convergence flag (0 = success, nonzero values are error codes) relative residual b − Ax / b number of iterations a vector listing the residual norm b − Ax at each iteration
The optional output variables are
The gmres method includes one additional input, [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,...) The method restarts every restart iterations, or if restart=[], the method does not restart.
19
Which method works best depends on the particular problem. This diagram (adapted from Demmel [3]) provides a reasonable starting point. A definite?
No Yes No Yes
Try pcg Try minres Try lsqr Try qmr
Yes No
A symmetric?
No
Is A wellconditioned? Is A wellYes conditioned?
Yes Yes No
A
available? No Is storage expensive?
Try cgs or bicgstab Try gmres
Functional Representation
Rather than storing the matrix explicitly, it is often convenient to represent A as a linear function acting on x. For example, 1 x1 1 1 x2 . . A*x = . .. . . . . . 1 1 · · · 1 xn can be expressed as A*x = cumsum(x). In addition to being conceptually convenient, functional representation avoids the memory overhead of storing the matrix. All the iterative solvers support functional representation, where a function handle afun is given in place of matrix A. The function afun(x) should accept a vector x and return A*x. The help for lsqr (“help lsqr”) has a functional example where A is the central differencing operator. The bicg, lsqr, and qmr methods also require the transpose: afun(x,’notransp’) should return A*x and afun(x,’transp’) should return A’*x. Similarly, the preconditioner may be given as a function handle mfun.
Special Matrices
For certain special matrix structures, the Ax = b problem can be solved very quickly.
Circulant matrices. Matrices corresponding to cyclic convolution Ax = h ∗ x are diagonalized in the Fourier domain, and can be solved by x = ifft(fft(b)./fft(h)); Triangular and banded. Triangular matrices and diagonally-dominant banded matrices are solved efficiently by sparse LU factorization: [L,U] = lu(sparse(A)); x = U\(L\b); Poisson problem. If A is a finite difference approximation of the Laplacian, the problem is efficiently solved by multigrid methods (e.g., web search for “matlab multigrid”).
20
8.2 Inlined PCG
As discussed in Section 6, it is sometimes advantageous to inline a simple piece of code to tailor it to a specific application. While most iterative methods are too lengthy for inlining, conjugate gradients has a short implementation. The value of inlining pcg is that calls to afun and mfun may be inlined. Beware that in the following code, A must be symmetric positive definite—otherwise, it will not work!
%% Use PCG to tol = 1e−6; maxit = 20; x = x0; solve Ax = b %% % Convergence tolerance % Maximum number of iterations % Set x to initial guess (if no guess is available, use x0 = 0)
Compute w = A x
r = b − w; % r = residual vector
Set p = r, or if using a preconditioner, p = M −1 r
d = r'*p; for iter = 1:maxit if norm(p) < tol break; end
% PCG converged successfully (small change in x)
Compute w = A p
alpha = e/(p'*w); x = x + alpha*p; r = r − alpha*w;
Set w = r, or if using a preconditioner, w = M −1 r
if norm(w) < tol && norm(r) < tol break; % PCG converged successfully (small residual) end dlast = d; d = r'*w; p = w + (d/dlast)*p; end
After the loop, the final solution is in x. The method converged successfully with accuracy tol if it encountered a break statement, otherwise, it failed or needed more iterations. The relative residual may be computed by relres = norm(r)/norm(b) after the loop and resvec by storing the value of norm(r) each iteration.
21
9
Numerical Integration
b
Quadrature formulas are numerical approximation of integrals of the form f (x) dx ≈
a k
wk f (xk ),
where the xk are called the nodes or abscissas and wk are the associated weights. Simpson’s rule is
b
f (x) dx ≈
a
h [f (a) + 4f (a + h) + f (b)] , 3
h=
b−a . 2
h 4h h 3, 3 , 3.
Simpson’s rule is a quadrature formula with nodes a, a + h, b and node weights Matlab has several functions for quadrature in one dimension: quad quadl quadgk quadv adaptive adaptive adaptive adaptive Simpson Gauss-Lobatto Gauss-Kronrod Simpson
Better for low accuracy and nonsmooth integrands Higher accuracy with smoother integrands Oscillatory integrands and high accuracy Vectorized for multi-valued integrands
The quad- functions are robust and precise, however, they are not very efficient. They use an adaptive refinement procedure, and while this reduces the number of function calls, they gain little from vectorization and incur significant overhead. If an application requires approximating an integral whose integrand can be efficiently vectorized, using nonadaptive quadrature may improve speed.
for n = −20:20 % Compute Fourier series coefficients c(n + 21) = quad(@(x)(exp(sin(x).ˆ6).*exp(−1i*x*n)),0,pi,1e−4); end
This code runs in 5.16 seconds. In place of quad, using Simpson’s composite rule with N = 199 nodes yields results with comparable accuracy and allows for vectorized computation. Since the integrals are all over the same interval, the nodes and weights only need to be constructed once.
N = 199; h = pi/(N−1); x = (0:h:pi).'; w = ones(1,N); w(2:2:N−1) = 4;
w(3:2:N−2) = 2;
w = w*h/3;
% Nodes % Weights
for n = −20:20 c(n + 21) = w * ( exp(sin(x).ˆ6).*exp(−1i*x*n) ); end
This version of the code runs in 0.02 seconds (200 times faster). The quadrature is performed by the dot product multiplication with w. It can be further optimized by replacing the for loop with one vector-matrix multiply:
[n,x] = meshgrid(−20:20, 0:h:pi); c = w * ( exp(sin(x).ˆ6).*exp(−1i*x.*n) );
22
For this example, quadv can be used on a multi-valued integrand with similar accuracy and speed,
n = −20:20; c = quadv(@(x)exp(sin(x).ˆ6).*exp(−1i.*x.*n),0,pi,1e−4);
9.1 One-Dimensional Integration
b a
f (x) dx is approximated by composite Simpson’s rule with
= = = = (b − a)/(N−1); (a:h:b).'; ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; w * f(x); % Approximately evaluate the integral
h x w I
where N is an odd integer specifying the number of nodes. A good higher-order choice is 4th-order Gauss-Lobatto [4] (as used by quadl), based on
1 −1 1 1 1 f (x) dx ≈ 6 f (−1) + 5 f (− √5 ) + 5 f ( √5 ) + 1 f (1). 6 6 6
N h d x w I
= = = = = =
max(3*round((N−1)/3),3) + 1; % Adjust N to the closest valid choice (b − a)/(N−1); (3/sqrt(5) − 1)*h/2; (a:h:b).'; x(2:3:N−2) = x(2:3:N−2) − d; x(3:3:N−1) = x(3:3:N−1) + d; ones(1,N); w(4:3:N−3) = 2; w([2:3:N−2,3:3:N−1]) = 5; w = w*h/4; w * f(x); % Approximately evaluate the integral
The number of nodes N must be such that (N − 1)/3 is an integer. If not, the first line adjusts N to the closest valid choice. It is usually more accurate than Simpson’s rule when f has six continuous derivatives, f ∈ C 6 (a, b). A disadvantage of this nonadaptive approach is that the accuracy of the result is only indirectly controlled by the parameter N. To guarantee a desired accuracy, either use a generously large value for N or, if possible, determine the error bounds [6, 7] Simpson’s rule error 4th-order Gauss-Lobatto error ≤ ≤ (b − a)h4 max f (4) (η) provided f ∈ C 4 (a, b), a≤η≤b 180 27(b − a)h6 max f (6) (η) provided f ∈ C 6 (a, b), a≤η≤b 56000
b−a where h = N −1 . Note that these bounds are valid only when the integrand is sufficiently differentiable: f must have four continuous derivatives for the Simpson’s rule error bound, and six continuous derivatives for 4th-order Gauss-Lobatto.
23
Composite Simpson’s rule is a fast and effective default choice. But depending on the situation, other methods may be faster and more accurate:
If the integrand is expensive to evaluate, or if high accuracy is required and the error bounds above are too difficult to compute, use one of Matlab’s adaptive methods. Otherwise, consider composite methods. Use higher-order methods (like Gauss-Lobatto/quadl) for very smooth integrands and lower-order methods for less smooth integrands. Use the substitution u =
1 1−x
or Gauss-Laguerre quadrature for infinite-domain integrals like
∞ . 0
9.2 Multidimensional Integration
An approach for evaluating double integrals of the form a c f (x, y) dy dx is to apply one-dimensional b quadrature to the outer integral a F (x) dx and then for each x use one-dimensional quadrature over d the inner dimension to approximate F (x) = c f (x, y) dy. The following code does this with composite Simpson’s rule with Nx×Ny nodes:
%%% Construct Simpson nodes and weights over x %%% h = (b − a)/(Nx−1); x = (a:h:b).'; wx = ones(1,Nx); wx(2:2:Nx−1) = 4; wx(3:2:Nx−2) = 2; %%% Construct Simpson nodes and weights over y %%% h = (d − c)/(Ny−1); y = (c:h:d).'; wy = ones(1,Ny); wy(2:2:Ny−1) = 4; wy(3:2:Ny−2) = 2; %%% Combine for two−dimensional integration %%% [x,y] = meshgrid(x,y); x = x(:); y = y(:); w = wy.'*wx; w = w(:).'; I = w * f(x,y); % Approximately evaluate the integral
b d
wx = w*h/3;
wy = w*h/3;
Similarly for three-dimensional integrals, the weights are combined with
[x,y,z] = meshgrid(x,y,z); w = wy.'*wx; w = w(:)*wz; x = x(:); y = y(:); w = w(:).'; z = z(:);
When the integration region is complicated or of high dimension, Monte Carlo integration techniques are appropriate. The disadvantage is that an N -point Monte Carlo quadrature has error on the order O( √1 ), so many points are necessary even for moderate accuracy. Suppose that N points, N x1 , x2 , . . . , xN , are uniformly randomly selected in a multidimensional region (or volume) Ω. Then f dV ≈
Ω Ω
dV N
N
f (xn ).
n=1
24
To integrate a complicated region W that is difficult to sample uniformly, find an easier region Ω that contains W and can be sampled [5]. Then f dV =
W Ω
f · χW dV ≈
Ω
dV N
N
f (xn )χW (xn ),
n=1
χW (x) =
1, 0,
x ∈ W, x ∈ W. /
χW (x) is the indicator function of W : χW (x) = 1 when x is within region W and χW (x) = 0 otherwise. Multiplying the integrand by χW sets contributions from outside of W to zero. For example, consider finding the center of mass of the shape W defined by cos 2
2 2
x2 + y 2 x ≤ y
and x + y ≤ 4. Given the integrals M = W dA, Mx = W x dA, and My = W y dA, the center of M mass is ( Mx , My ). The region is contained in the rectangle Ω defined by −2 ≤ x ≤ 2 and −2 ≤ y ≤ 2. M The following code estimates M , Mx , and My :
%%% Uniformly randomly sample points (x,y) in Ω %%% x = 4*rand(N,1) − 2; y = 4*rand(N,1) − 2; %%% Restrict the points to region W %%% i = (cos(2*sqrt(x.ˆ2 + y.ˆ2)).*x <= y & x.ˆ2 + y.ˆ2 <= 4); x = x(i); y = y(i); %%% Approximately evaluate the integrals %%% area = 4*4; % The area of rectangle Ω M = (area/N) * length(x); Mx = (area/N) * sum(x); My = (area/N) * sum(y);
2
0
-2 -2
0
2
Region W sampled with N = 1500. The center of mass is ≈ (0.5, 0.7).
More generally, if W is a two-dimensional region contained in the rectangle defined by a ≤ x ≤ b and c ≤ y ≤ d, the following code approximates W f dA:
x y i x = = = = a + (b−a)*rand(N,1); c + (d−c)*rand(N,1); logical(indicatorW(x,y)); x(i); y = y(i);
area = (b−a)*(d−c); I = (area/N) * sum(f(x,y));
% Approximately evaluate the integral
where indicatorW(x,y) is the indicator function χW (x, y) for region W . For refinements and variations of Monte Carlo integration, see for example [1].
25
10
Signal Processing
Even without the Signal Processing Toolbox, Matlab is quite capable in signal processing computations. This section lists code snippets to perform several common operations efficiently.
Moving-average filter
To compute an N-sample moving average of x with zero padding:
y = filter(ones(N,1)/N,1,x);
For large N, it is faster to use
y = cumsum(x)/N; y(N+1:end) = y(N+1:end) − y(1:end−N);
Locating zero-crossings and extrema
To obtain the indices where signal x crosses zero:
i = find(diff(sign(x))); % The kth zero−crossing lies between x(i(k)) and x(i(k)+1)
Linear interpolation can be used for subsample estimates of zero-crossing locations:
i = find(diff(sign(x))); i = i − x(i)./(x(i+1) − x(i)); % Linear interpolation
Since local maximum and minimum points of a signal have zero derivative, their locations can be estimated from the zero-crossings of diff(x), provided the signal is sampled with sufficiently fine resolution. For a coarsely sampled signal, a better estimate is
iMax = find(sign(x(2:end−1)−x(1:end−2)) + sign(x(2:end−1)−x(3:end)) > iMin = find(sign(x(2:end−1)−x(1:end−2)) + sign(x(2:end−1)−x(3:end)) < ... 0) + 1; ... 0) + 1;
FFT-based convolution
This line performs FFT-based circular convolution, equivalent to y = filter(b,1,x) except near the boundaries, provided that length(b) < length(x):
y = ifft(fft(b,length(x)).*fft(x));
For FFT-based zero-padded convolution, equivalent to y = filter(b,1,x),
N = length(x)+length(b)−1; y = ifft(fft(b,N).*fft(x,N)); y = y(1:length(x));
In both code snippets above, y will be complex-valued, even if x and b are both real. To force y to be real, follow the computation with y = real(y). If you have the Signal Processing Toolbox, it is faster to use fftfilt for FFT-based, zero-padded filtering. 26
Noncausal filtering and other boundary extensions
For its intended use, the filter command is limited to causal filters, that is, filters that do not involve “future” values to the right of the current sample. Furthermore, filter is limited to zero-padded boundary extension; data beyond the array boundaries are assumed zero as if the data were padded with zeros. For two-tap filters, noncausal filtering and other boundary extensions are possible through filter’s fourth initial condition argument. Given a boundary extension value padLeft for x(0), the filter yn = b1 xn + b2 xn−1 may be implemented as
y = filter(b,1,x,padLeft*b(2));
Similarly, given a boundary extension value padRight for x(end+1), the filter yn = b1 xn+1 + b2 xn is implemented as
y = filter(b,1,[x(2:end),padRight],x(1)*b(2));
Choices for padLeft and padRight for various boundary extensions are Boundary extension Periodic Whole-sample symmetric Half-sample symmetric Antisymmetric padLeft x(end) x(2) x(1) 2*x(1)-x(2) padRight x(1) x(end-1) x(end) 2*x(end)-x(end-1)
It is in principle possible to use a similar approach for longer filters, but ironically, computing the initial condition itself requires filtering. To implement noncausal filtering and filtering with other boundary handling, it is usually fastest to pad the input signal, apply filter, and then truncate the result.
Upsampling and Downsampling
Upsample x (insert zeros) by factor p:
y = zeros(length(x)*p−p+1,1); y(1:p:length(x)*p) = x; % For trailing zeros, use y = zeros(length(x)*p,1);
Downsample x by factor p, where 1 ≤ q ≤ p:
y = x(q:p:end);
27
Haar Wavelet
This code performs K stages of the Haar wavelet transform on the input data x to produce wavelet coefficients y. The input array x must have length divisible by 2K . Forward transform:
y = x(:); N = length(y); for k = 1:K tmp = y(1:2:N) + y(2:2:N); y([1:N/2,N/2+1:N]) = ... [tmp;y(1:2:N) − 0.5*tmp]/sqrt(2); N = N/2; end
Inverse transform:
x = y(:); N = length(x)*pow2(−K); for k = 1:K N = N*2; tmp = x(N/2+1:N) + 0.5*x(1:N/2); x([1:2:N,2:2:N]) = ... [tmp;x(1:N/2) − tmp]*sqrt(2); end
Daubechies 4-Tap Wavelet
This code implements the Daubechies 4-tap wavelet in lifting scheme form [2]. The input array x must have length divisible by 2K . Filtering is done with symmetric boundary handling. Forward transform:
U = 0.25*[2−sqrt(3),−sqrt(3)]; ScaleS = (sqrt(3) − 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2); y = x(:); N = length(y); for k = 1:K N = N/2; y1 = y(1:2:2*N); y2 = y(2:2:2*N) + sqrt(3)*y1; y1 = y1 + filter(U,1,... y2([2:N,max(N−1,1)]),y2(1)*U(2)); y(1:2*N) = [ScaleS*... (y2 − y1([min(2,N),1:N−1]));ScaleD*y1]; end
Inverse transform:
U = 0.25*[2−sqrt(3),−sqrt(3)]; ScaleS = (sqrt(3) − 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2); x = y(:); N = length(x)*pow2(−K); for k y1 y2 y1 = 1:K = x(N+1:2*N)/ScaleD; = x(1:N)/ScaleS + y1([min(2,N),1:N−1]); = y1 − filter(U,1,... y2([2:N,max(N−1,1)]),y2(1)*U(2)); x([1:2:2*N,2:2:2*N]) = [y1;y2 − sqrt(3)*y1]; N = 2*N; end
To use periodic boundary handling rather than symmetric boundary handling, change appearances of [2:N,max(N-1,1)] to [2:N,1] and [min(2,N),1:N-1] to [N,1:N-1].
28
11
Miscellaneous
Clip a value without using if statements
To clip (or clamp, bound, or saturate) a value to within a range, the straightforward implementation is
if x < x = elseif x = end lowerBound lowerBound; x > upperBound upperBound;
Unfortunately, this is slow. A faster method is to use the min and max functions:
x = max(x,lowerBound); x = min(x,upperBound); % Clip elements from below, force x >= lowerBound % Clip elements from above, force x <= upperBound
Moreover, it works per-element if x a matrix of any size.
Convert any array into a column vector
It is often useful to force an array to be a column vector, for example, when writing a function expecting a column vector as an input. This simple trick will convert the input array of any shape and number of dimensions to a column vector, even if the input is already a column vector.
x = x(:); % Convert x to a column vector
Similarly, x = x(:,:) reduces an N-d array to 2D, where the number of rows stays the same.
Find the min/max of a matrix or N-d array
Given a matrix input (with no other inputs), the min and max functions operate along the columns, finding the extreme element in each column. To find the extreme element over the entire matrix, one way for a two-dimensional matrix is min(min(A)). A better method that works regardless of the number of dimensions and determines the extreme element’s location is
% Find the minimum element in A % The minimum value is MinValue, the index is MinIndex MinSub = ind2sub(size(A), MinIndex); % Convert MinIndex to subscripts [MinValue, MinIndex] = min(A(:));
The minimum element is A(MinIndex) or A(MinSub(1), MinSub(2), . . .) as a subscript reference. (Similarly, replace min with max for maximum value.)
29
Flood filling
Flood filling, like the “bucket” tool in image editors, can be elegantly written as a recursive function:
function I = flood1(I,c,x,y) % Flood fills image I from point (x,y) with color c. c2 = I(y,x); I(y,x) = c; if if if if x x y y > < > < 1 size(I,2) 1 size(I,1) & & & & I(y,x−1) I(y,x+1) I(y−1,x) I(y+1,x) == == == == c2 c2 c2 c2 & & & & I(y,x−1) I(y,x+1) I(y−1,x) I(y+1,x) ˜= ˜= ˜= ˜= c, c, c, c, I I I I = = = = flood1(I,c,x−1,y); flood1(I,c,x+1,y); flood1(I,c,x,y−1); flood1(I,c,x,y+1); end end end end
Being a highly recursive function, this is inefficient in Matlab. The following code is faster:
function I = flood2(I,c,x,y) % Flood fills image I from point (x,y) with color c. LastFlood = zeros(size(I)); Flood = LastFlood; Flood(y,x) = 1; Mask = (I == I(y,x)); FloodFilter = [0,1,0; 1,1,1; 0,1,0]; while any(LastFlood(:) ˜= Flood(:)) LastFlood = Flood; Flood = conv2(Flood,FloodFilter,'same') & Mask; end I(find(Flood)) = c;
The key is the conv2 two-dimensional convolution function. Flood filling a 40×40-pixel region takes 1.168 seconds with flood1 and 0.067 seconds with flood2.
Vectorized use of set on GUI objects
While not a speed improvement, a trick for cleaner code is vectorized use of set on GUI objects. A serious graphical user interface (GUI) can have dozens of objects and many properties to configure. For example, to define three edit boxes with white text background and left text alignment:
uicontrol('Units', 'normalized', 'Position', [0.1,0.9,0.7,0.05], ... 'HorizontalAlignment', 'left', 'Style', 'edit', 'BackgroundColor', [1,1,1]); uicontrol('Units', 'normalized', 'Position', [0.1,0.8,0.7,0.05], ... 'HorizontalAlignment', 'left', 'Style', 'edit', 'BackgroundColor', [1,1,1]); uicontrol('Units', 'normalized', 'Position', [0.1,0.7,0.7,0.05], ... 'HorizontalAlignment', 'left', 'Style', 'edit', 'BackgroundColor', [1,1,1]);
30
A vectorized call to set reduces this to
h(1) = uicontrol('Units', 'normalized', 'Position', [0.1,0.9,0.7,0.05]); h(2) = uicontrol('Units', 'normalized', 'Position', [0.1,0.8,0.7,0.05]); h(3) = uicontrol('Units', 'normalized', 'Position', [0.1,0.7,0.7,0.05]); set(h, 'HorizontalAlignment', 'left', 'Style', 'edit','BackgroundColor', [1,1,1]);
12
Further Reading
For more reading on vectorization, see the MathWorks vectorization guide: http://www.mathworks.com/support/tech-notes/1100/1109.shtml For further details and examples on JIT Acceleration, see for example http://www.mathworks.com/access/helpdesk_r13/help/techdoc/matlab_prog/ch7_perf.html For a thorough guide on efficient array manipulation, see MATLAB array manipulation tips and tricks: http://home.online.no/~pjacklam/matlab/doc/mtt For numerical methods with Matlab, see Numerical Computing with MATLAB : http://www.mathworks.com/moler
MEX
In a coding situation that cannot be optimized any further, keep in mind that the Matlab language is intended primarily for easy prototyping rather than high-speed computation. In some cases, an appropriate solution is Matlab Executable (MEX) external interface functions. With a C or Fortran compiler, it is possible to produce a MEX function that can be called from within Matlab in the same as an M-function. The typical speed improvement over equivalent M-code is easily ten-fold. Installations of Matlab include an example of MEX under directory <MATLAB>/extern/examples/mex. In this directory, a function “yprime” is written as M-code (yprime.m) and equivalent C (yprime.c) and Fortran (yprime.f) MEX source code. For more information, see the MathWorks MEX-files guide http://www.mathworks.com/support/tech-notes/1600/1605.html For information on compiling MEX files with the GNU compiler collection (GCC), see http://gnumex.sourceforge.net Beware, MEX is an ambitious alternative to M-code. Writing MEX files requires learning some low-level details about Matlab, and takes more time to develop than M-code.
31
Matlab Compiler
The Matlab Compiler generates a standalone executable from M-code. For product information and documentation, see http://www.mathworks.com/access/helpdesk/help/toolbox/compiler
References
[1] A. Bielajew. “The Fundamentals of the Monte Carlo Method for Neutral and Charged Particle Transport.” University of Michigan, class notes. Available online at http://www-personal.engin.umich.edu/~bielajew/MCBook/book.pdf [2] I. Daubechies and W. Sweldens. “Factoring Wavelet Transforms into Lifting Steps.” Sep. 1996. [3] J. W. Demmel. Applied Numerical Linear Algebra. SIAM, 1997. [4] W. Gander and W. Gautschi. “Adaptive Quadrature–Revisited,” BIT, vol. 40, pp. 84-101, 2000. [5] W. Press, B. Flannery, S. Teukolsky and W. Vetterling. Numerical Recipes. Cambridge University Press, 1986. [6] E. Weisstein. “Lobatto Quadrature.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/LobattoQuadrature.html [7] E. Weisstein. “Simpson’s Rule.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/SimpsonsRule.html
32
Writing fast Matlab code (ang)
Pascal Getreuer, February 2009
Contents
Page
1 Introduction 2 The Profiler 3 Array Preallocation 4 JIT Acceleration 5 Vectorization 6 Inlining Simple Functions 7 Referencing Operations 8 Solving Ax = b 9 Numerical Integration 10 Signal Processing 11 Miscellaneous 12 Further Reading
1 3 5 7 9 14 16 19 22 26 29 31
1
Introduction
Matlab is a prototyping environment, meaning it focuses on the ease of development with language flexibility, interactive debugging, and other conveniences lacking in performance-oriented languages like C and Fortran. While Matlab may not be as fast as C, there are ways to bring it closer.
Disclaimer This is not a beginner’s tutorial to Matlab, but a tutorial on
performance. The methods discussed here are generally fast, but no claim is made on what is fastest. Use at your own risk.
Before beginning our main topic, let’s step back for a moment and ask ourselves what we really want. Do we want the fastest possible code?—if so, we should switch to a performance language like C. Probably more accurately, we want to spend less time total from developing, debugging, running, and until finally obtaining results. This section gives some tips on working smart, not hard.
M-Lint Messages
In more recent versions of Matlab, the integrated editor automatically gives feedback on possible code errors and suggested optimizations. These messages are helpful in improving a program and in learning more about Matlab.
for k = 1:NumTrials r = rand; x(k) = runsim(r); x might be growing inside a loop. Consider preallocating for speed. end hist(x); % Plot histogram of simulation outputs
Hold the mouse cursor over underlined code to read a message. Or, see all M-Lint messages at once ::::::::::::::: with Tools → M-Lint → Show M-Lint Report.
Organization
Use a separate folder for each project.
A separate folder for each project keeps related things together and simplifies making copies and backups of project files.
Write header comments, especially H1.
The first line of the header comment is called the H1 comment. Type help(cd) to get a listing of all project files and their H1 comments.
1
Save frequent console commands as a script.
If you find yourself repeating certain commands on the console, save them as a script. Less typing means fewer opportunities for typos.
Avoid losing data
Don’t use clear all in a script.
This is an unfortunate common practice—any important variables in the base workspace will be irretrievably lost.
Beware of clobber.
“File clobber” refers to the kind of data loss when a file is accidentally overwritten with another one having the same filename. This phenomenon can occur with variables as well:
>> result = big operation(input1); ... >> result = big operation(input2); % Store the result of a time−consuming operation % Run again with another input
Variable result was clobbered and the first output was lost.
Beware of what can crash MATLAB.
While Matlab is generally reliable, crashes are possible when using third-party MEX functions or extremely memory-intensive operations, for example, with video and very large arrays. Now with good working habits covered, we begin our discussion of writing fast Matlab code. The rest of this article is organized by topic, first on techniques that are useful in general application, and next on specific computational topics (table of contents is on the first page).
2
2
The Profiler
Matlab 5.0 (R10) and newer versions include a tool called the profiler that helps identify bottlenecks in a program. Use it with the profile command: profile profile profile profile on off clear report Turn the profiler on Turn it back off Clear profile statistics View the results from the profiler
For example, consider profiling 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 and clear the profiler, run the function, and then view the profile report:
>> profile on, profile clear >> example1(5000); >> profile report
There is a slight parsing overhead when running code for the first time; run the test code twice and time the second run. The profiler report command shows a report. Depending on the system, profiler results may differ from this example.
MATLAB Profile Report: Summary
Report generated 30-Jul-2004 16:57:01 Total recorded time: Number of M-functions: Clock precision: 3.09 s 4 0.016 s
3
Function List
Name
example1 gammaln profile profreport
Time
3.09 0.73 0.00 0.00 100.0% 23.7% 0.0% 0.0%
Time
1 3562 1 1
Time/call
3.094000 0.000206 0.000000 0.000000
Self time
2.36 0.73 0.00 0.00 76.3% 23.7% 0.0% 0.0%
Location
~/example1.m ../toolbox/matlab/specfun/gammaln.m ../toolbox/matlab/general/profile.m ../toolbox/matlab/general/profreport.m
Clicking the “example1” link gives more details: Lines where the most time was spent
Line Number Code result(k) = sin(k/50); result(k) = gammaln(k); if result(k) < -0.9 Calls Total Time % Time
4 7 6
Totals
5000 721 5000
2.11 0.84 0.14 3.09
s s s s
68% 27% 5% 100%
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. Another helpful section of the profile report is “M-Lint Results,” which gives feedback from the M-Lint code analyzer. Possible errors and suggestions are listed here. M-Lint results
Line number Message ‘result’ might be growing inside a loop. Consider preallocating for speed. ‘result’ might be growing inside a loop. Consider preallocating for speed.
4 7
(Preallocation is discussed in the next section.) The profiler has limited time resolution, so to profile a piece of code that runs too quickly, run the test code multiple times with a loop. Adjust the number of loop iterations so that the time it takes to run is noticeable. More iterations yields better time resolution. The profiler is an essential tool for identifying bottlenecks and per-statement analysis, however, for more accurate timing of a piece of code, use the tic/toc stopwatch timer.
>> tic; example1(5000); toc; Elapsed time is 3.082055 seconds
For serious benchmarking, also close your web browser, anti-virus, and other background processes that may be taking CPU cycles.
4
3
>> a = 2 a = 2
Array Preallocation
Matlab’s matrix variables have the ability to dynamically augment rows and columns. For example,
>> 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 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
This code takes 0.47 seconds to run. 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.
a = zeros(1,8000); b = zeros(1,8000); a(1) = 1; b(1) = 0; % Preallocation
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.
5
What if the final array size can vary? Here is an example:
a = zeros(1,10000); count = 0; for k = 1:10000 v = exp(rand*rand); if v > 0.5 % Conditionally add to array count = count + 1; a(count) = v; end end a = a(1:count); % Trim extra zeros from the results % Preallocate
The average run time of this program is 0.42 seconds without preallocation and 0.18 seconds with it. Preallocation is also beneficial for cell arrays, using the cell command to create a cell array of the desired size.
6
4
JIT Acceleration
Matlab 6.5 (R13) and later feature the Just-In-Time (JIT) Accelerator for improving the speed of M-functions, particularly with loops. By knowing a few things about the accelerator, you can improve its performance. The JIT Accelerator is enabled by default. To disable it, type “feature accel off” in the console, and “feature accel on” to enable it again. As of Matlab R2008b, only a subset of the Matlab language is supported for acceleration. Upon encountering an unsupported feature, acceleration processing falls back to non-accelerated evaluation. Acceleration is most effective when significant contiguous portions of code are supported.
Data types: Code must use supported data types for acceleration: double (both real and complex), logical, char, int8–32, uint8–32. Some struct, cell, classdef, and function handle usage is supported. Sparse arrays are not accelerated. Array shapes: Array shapes of any size with 3 or fewer dimensions are supported. Changing the shape or data type of an array interrupts acceleration. A few limited situations with 4D arrays are accelerated. Function calls: Calls to built-in functions and M-functions are accelerated. Calling MEX functions and Java interrupts acceleration. (See also page 14 on inlining simple functions.) Conditionals and loops: The conditional statements if, elseif, and simple switch statements are supported if the conditional expression evaluates to a scalar. Loops of the form for k=a:b, for k=a:b:c, and while loops are accelerated if all code within the loop is supported.
In-place computation
Introduced in Matlab 7.3 (R2006b), the element-wise operators (+, .*, etc.) and some other functions can be computed in-place. That is, a computation like
x = 5*sqrt(x.ˆ2 + 1);
is handled internally without needing temporary storage for accumulating the result. An M-function can also be computed in-place if its output argument matches one of the input arguments.
x = myfun(x); function x = myfun(x) x = 5*sqrt(x.ˆ2 + 1); return;
To enable in-place computation, the in-place operation must be within an M-function (and for an inplace function, the function itself must be called within an M-function). Currently, there is no support for in-place computation with MEX-functions.
7
Multithreaded Computation
Matlab 7.4 (R2007a) introduced multithreaded computation for multicore and multiprocessor computers. Multithreaded computation accelerates some per-element functions when applied to large arrays (for example, .^, sin, exp) and certain linear algebra functions in the BLAS library. To enable it, select File → Preferences → General → Multithreading and select “Enable multithreaded computation.” Further control over parallel computation is possible with the Parallel Computing Toolbox using parfor and spmd.
JIT-Accelerated Example
For example, the following loop-heavy code is supported for acceleration:
function B = bilateral(A,sd,sr,R) % The bilateral image denoising filter B = zeros(size(A)); for i = 1:size(A,1) for j = 1:size(A,2) zsum = 0; for m = −R:R if i+m >= 1 && i+m <= size(A,1) for n = −R:R if j+n >= 1 && j+n <= size(A,2) z = exp(−(A(i+m,j+n) − A(i,j))ˆ2/(2*sdˆ2))... * exp(−(mˆ2 + nˆ2)/(2*srˆ2)); zsum = zsum + z; B(i,j) = B(i,j) + z*A(i+m,j+n); end end end end B(i,j) = B(i,j)/zsum; end end
For a 128 × 128 input image and R = 3, the run time is 53.3 seconds without acceleration and 0.68 seconds with acceleration.
8
5
Vectorization
A computation is vectorized by taking advantage of vector operations. A variety of programming situations can be vectorized, and often improving speed to 10 times faster or even better. Vectorization is one of the most general and effective techniques for writing fast M-code.
5.1 Vectorized Computations
Many standard Matlab 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 3 >> abs([0,1,2,−5,−6,−7]) ans = 2 4 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, its distance from the origin is computed and stored in d. For speed, array d is preallocated (see Section 3). 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); d = min(d); % Compute distance for every point % 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. (A further improvement: it is equivalent to compute d = sqrt(min(x.^2 + y.^2 + z.^2))).
9
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, accumarray
5.2 Vectorized Logic
The previous section shows how to vectorize pure computation. But bottleneck code often involves conditional 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 values. How is this useful? Matlab has a few powerful functions for operating on logical arrays:
>> find([1,5,3] < [2,2,4]) % Find indices of nonzero elements ans = 1
3 % True if any element of a vector is nonzero % (or per−column for a matrix)
>> any([1,5,3] < [2,2,4]) ans = 1 >> all([1,5,3] < [2,2,4]) ans = 0
% True if all elements of a vector are nonzero % (or per−column for a matrix)
Vectorized logic also works on arrays.
>> find(eye(3) == 1) ans = 1 5 9 % (element locations given as indices) % Find which elements == 1 in the 3x3 identity matrix
By default, find returns element locations as indices (see page 16 for indices vs. subscripts).
10
Example 1: Vector Normalization
To normalize a single vector v to unit length, one can use v = v/norm(v). However, to use norm for set of vectors v(:,1), v(:,2), . . . requires computing v(:,k)/norm(v(:,k)) in a loop. Alternatively,
vMag = sqrt(sum(v.ˆ2)); v = v./vMag(ones(1,size(v,1)),:);
Example 2: 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)); x(i) = []; % Find bad elements in x % and delete them
Alternatively,
i = find(˜isnan(x) & ˜isinf(x)); x = x(i); % Find elements that are not NaN and not infinite % 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 3: Piecewise functions
The sinc function has a piecewise definition, sinc(x) = sin(x)/x, 1, x=0 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)); i = find(x ˜= 0); y(i) = sin(x(i)) ./ x(i); % Set y to all ones, sinc(0) = 1 % Find nonzero x values % Compute sinc where x ˜= 0
A concise alternative is y = (sin(x) + (x == 0))./(x + (x == 0)).
11
Example 4: 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 1 1 2 2 2 3 3 3 4 4 4 5 5 5 y = 1 2 3 1 2 3 1 2 3 1 2 3 1 2 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 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); % Spiral centered at (75,50) Img = sin(r/3 + th); imagesc(Img); colormap(hot); axis image, axis off
% Convert to polar
12
Example 5: Polynomial interpolation
Given n points x1 , . . . xn and corresponding function values y1 , . . . yn , the coefficients c0 , . . . , cn−1 of the interpolating polynomial P (x) = cn−1 xn−1 + · · · + c1 x + c0 can be found by solving x1 n−1 x2 n−1 . . . n−1 xn
x1 n−2 x2 n−2 xn n−2
··· ··· ···
x1 2 x2 2 xn 2
x1 x2 xn
1 1 . . . 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(:); y = y(:); n = length(x); % Make sure x and y are both column vectors % n = Number of points
% Construct the matrix on the left−hand−side 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×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(2) x(2) · · · . . . x(n) x(n) · · · x(1) x(2) . . . x(n) powMatrix = n−1 n−2 n−1 n−2 . . . n−1 n−2 ··· ··· ··· 0 0 . . . 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). The code above is only an example—use the standard polyfit function for serious use.
13
6
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 determines 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 is 0.047 seconds, nearly 9 times faster. Furthermore, by inlining median, it can be specifically tailored to evaluating 5-sample medians. But this is only an example; if the Signal Processing Toolbox is available, y = medfilt1(x,5) is much faster. A surprising number of Matlab’s functions are implemented as M-files, of which many can be inlined in just a few lines. In Matlab 5.3 (R11), the following functions are worth inlining: conv, cross, fft2, fliplr, flipud, ifft, ifft2, ifftn, ind2sub, ismember, linspace, logspace, mean, median, meshgrid, poly, polyval, repmat, roots, rot90, setdiff, setxor, sortrows, std, sub2ind, union, unique, var The list above is for Matlab R11; some of these functions have since become built-in. Use the which command or try “edit function name” to determine whether a function is implemented as an M-file.
14
For example, in Matlab R11 and earlier, ifft was implemented by simply calling fft and conj. If x is a one-dimensional array, y = ifft(x) can be inlined with y = conj(fft(conj(x)))/length(x). Another example: b = unique(a) can be inlined with
b = sort(a(:)); b( b((1:end−1)') == b((2:end)') ) = [];
While repmat has the generality to operate on matrices, it is often only necessary to tile a vector or just a scalar. To repeat a column vector y over the columns n times,
A = y(:,ones(1,n)); % Equivalent to A = repmat(y,1,n);
To repeat a row vector x over the rows m times,
A = x(ones(1,m),:); % Equivalent to A = repmat(x,m,1);
To repeat a scalar s into an m×n matrix,
A = s(ones(m,n)); % Equivalent to A = repmat(s,m,n);
This method avoids the overhead of calling an M-file function. It is never slower than repmat (critics should note that repmat.m itself uses this method to construct mind and nind). For constructing matrices with constant value, there are other efficient methods, for example, s+zeros(m,n). Don’t go overboard. Inlining functions is only beneficial when the function is simple and when it is called frequently. Doing it unnecessarily obfuscates the code.
15
7
Referencing Operations
Referencing in Matlab is varied and powerful enough to deserve a section of discussion. Good understanding of referencing enables vectorizing a broader range of programming situations.
7.1 Subscripts vs. Indices
Subscripts are the most common method used to refer to matrix elements, for example, A(3,9) refers to row 3, column 9. Indices are an alternative referencing method. Consider a 10 × 10 matrix A. Internally, Matlab stores the matrix data linearly as a one-dimensional, 100-element array of data in column-major order. 1 2 3 . . . 10 11 21 12 22 13 23 . . . . . . 20 30 ··· 81 82 83 . . . 90 91 92 93 . . . 100
···
An index refers to an element’s position in this one-dimensional array. Using an index, A(83) also refers to the element on row 3, column 9. It is faster to scan down columns than over rows. Column-major order means that elements along a column are sequential in memory while elements along a row are further apart. Scanning down columns promotes cache efficiency. Conversion between subscripts and indices can be done with the sub2ind and ind2sub functions. However, because these are M-file functions rather than fast built-in operations, it is more efficient to compute conversions directly (also see Section 6). For a two-dimensional matrix A of size M×N, the conversions between subscript (i,j) and index (index) are A(i,j) ↔ A(i + (j-1)*M) A(index) ↔ A(rem(index-1,M)+1, floor((index-1)/M)+1) Indexing extends to N-D matrices as well, with indices increasing first through the columns, then through the rows, through the third dimension, and so on. Subscript notation extends intuitively, A(. . . , dim4, dim3, row, col).
7.2 Vectorized Subscripts
It is useful to work with submatrices rather than individual elements. This is done with a vector of indices or subscripts. If A is a two-dimensional matrix, a vector subscript reference has the syntax A(rowv, colv)
16
where rowv and colv are vectors. Both may be of any length and their elements may be in any order. If either is a matrix, it is reshaped to a vector. There is no difference between using row vectors or column vectors in vector subscripts. Let M = length(rowv) and N = length(colv). side∗ returns a submatrix of size M×N: A(rowv(1), colv(1)) A(rowv(1), colv(2)) A(rowv(2), colv(1)) A(rowv(2), colv(2)) . . . A(rowv(M), colv(1)) A(rowv(M), colv(2)) Then a vector subscripted matrix on the right-hand
A(rowv(1), colv(3)) · · · A(rowv(2), colv(3)) · · · A(rowv(M), colv(3)) · · ·
A(rowv(1), colv(N)) A(rowv(2), colv(N)) . . . A(rowv(M), colv(N))
A vector subscripted matrix can also appear on the left-hand side as A(rowv, colv) = scalar or M×N matrix If any elements in the destination reference are repeated, the source element with the greater index dominates. For example,
>> A = zeros(2); A([1,2],[2,2]) = [1,2;3,4] A = 0 0 2 4
Vector subscript references extend intuitively in higher dimensions.
7.3 Vector Indices
Multiple elements can also be referenced with vector indices. A(indexv) where indexv is an array of indices. On the right-hand side, a vector index reference returns a matrix the same size and shape as indexv. If indexv is a 3 × 4 matrix, then A(indexv) is the 3 × 4 matrix A(indexv(1, 1)) A(indexv(1, 2)) A(indexv(1, 3)) A(indexv(1, 4)) A(indexv) = A(indexv(2, 1)) A(indexv(2, 2)) A(indexv(2, 3)) A(indexv(2, 4)) A(indexv(3, 1)) A(indexv(3, 2)) A(indexv(3, 3)) A(indexv(3, 4)) While vector subscripts are limited to referring to block-shaped submatrices, vector indices can refer to any subset of elements. If a vector index reference is on the left-hand side, the right-hand side must return a matrix of the same size as indexv or a scalar. As with vector subscripts, ambiguous duplicate assignments use the value with greater source index.
∗ The
left- and right-hand sides (LHS and RHS) refer to the two sides of an assignment, “LHS = RHS.”
17
7.4 Reference Wildcards
Using the wildcard, :, in a subscript refers to an entire row or column. For example, A(:,1) refers to every row in column one–the entire first column. This can be combined with vector subscripts, A([2,4],:) refers to the second and fourth rows. When the wildcard is used in a vector index, the entire matrix is referenced. On the right-hand side, this always returns a column vector. A(:) = column vector This is frequently useful: for example, if a function input must be a row vector, the user’s input can be quickly reshaped into row vector with A(:).’ (make a column vector and transpose to a row vector). A(:) on the left-hand side assigns all the elements of A, but does not change its size. For example, A(:) = 8 changes all elements of matrix A to 8.
7.5 Logical Indexing
Given a logical array mask with the same size as A,
A(mask)
refers to all elements of A where the corresponding element of mask is 1 (true). It is equivalent to A(find(mask)). A common usage of logical indexing is to refer to elements based on a per-element condition, for example,
A(abs(A) < 1e−3) = 0
sets all elements with magnitude less than 10−3 to zero. Logical indexing is also useful to select nonrectangular sets of elements, for example,
A(logical(eye(size(A))))
references the diagonal elements of A. Note that for a right-hand side reference, it is faster to use diag(A) to get the diagonal of A.
7.6 Deleting Submatrices with [ ]
Elements in a matrix can be deleted by assigning the empty matrix. For example, A([3,5]) = [ ] deletes the third and fifth element from A. If this is done with index references, the matrix is reshaped into a row vector. It is also possible to delete with subscripts if all but one subscript are the wildcard. A(2,:) deletes the second row. Deletions like A(2,1) = [ ] or even A(2,1:end) = [ ] are illegal. = []
18
8
Solving Ax = b
Many problems involve solving a linear system, that is, solving the Ax = b problem b1 x1 a1,1 a1,2 · · · a1,n a2,1 a2,2 · · · a2,n x2 b2 . . . . = . . .. . . . . . . . . . . . bn xn an,1 an,2 · · · an,n The condition number cond(A) describes how sensitive the problem is to small perturbations. Matlab can efficiently approximate cond(A) with condest. Smaller condition number is better, and in this case A is said to be well-conditioned. An infinite condition number implies that A is singular—the problem cannot be solved, though it may be solved approximately with the pseudoinverse (x = pinv(A)*b). Small- to moderately-sized problems are solved efficiently by the backslash operator:
x = A\b; % Solves A*x = b
We shall focus on large problems.
8.1 Iterative Methods
For large problems, Matlab is well-equipped with iterative solvers: bicg, bicgstab, cgs, gmres, lsqr, minres, pcg, symmlq, qmr. The general calling syntax (except for gmres) is [x,flag,relres,iter,resvec] = method (A,b,tol,maxit,M1,M2,x0) where the required arguments are the matrix A, the right-hand side b, and the solution x. The minres, pcg, and symmlq methods require A to be symmetric (A = A’) and pcg additionally requires positive definiteness (all eigenvalues are positive). The optional input variables are tol maxit M1, M2 x0 flag relres iter resvec accuracy tolerance (default 10−6 ) maximum number of iterations (default 20) the method uses preconditioner M = M1*M2; it solves M −1 Ax = M −1 b initial guess convergence flag (0 = success, nonzero values are error codes) relative residual b − Ax / b number of iterations a vector listing the residual norm b − Ax at each iteration
The optional output variables are
The gmres method includes one additional input, [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,...) The method restarts every restart iterations, or if restart=[], the method does not restart.
19
Which method works best depends on the particular problem. This diagram (adapted from Demmel [3]) provides a reasonable starting point. A definite?
No Yes No Yes
Try pcg Try minres Try lsqr Try qmr
Yes No
A symmetric?
No
Is A wellconditioned? Is A wellYes conditioned?
Yes Yes No
A
available? No Is storage expensive?
Try cgs or bicgstab Try gmres
Functional Representation
Rather than storing the matrix explicitly, it is often convenient to represent A as a linear function acting on x. For example, 1 x1 1 1 x2 . . A*x = . .. . . . . . 1 1 · · · 1 xn can be expressed as A*x = cumsum(x). In addition to being conceptually convenient, functional representation avoids the memory overhead of storing the matrix. All the iterative solvers support functional representation, where a function handle afun is given in place of matrix A. The function afun(x) should accept a vector x and return A*x. The help for lsqr (“help lsqr”) has a functional example where A is the central differencing operator. The bicg, lsqr, and qmr methods also require the transpose: afun(x,’notransp’) should return A*x and afun(x,’transp’) should return A’*x. Similarly, the preconditioner may be given as a function handle mfun.
Special Matrices
For certain special matrix structures, the Ax = b problem can be solved very quickly.
Circulant matrices. Matrices corresponding to cyclic convolution Ax = h ∗ x are diagonalized in the Fourier domain, and can be solved by x = ifft(fft(b)./fft(h)); Triangular and banded. Triangular matrices and diagonally-dominant banded matrices are solved efficiently by sparse LU factorization: [L,U] = lu(sparse(A)); x = U\(L\b); Poisson problem. If A is a finite difference approximation of the Laplacian, the problem is efficiently solved by multigrid methods (e.g., web search for “matlab multigrid”).
20
8.2 Inlined PCG
As discussed in Section 6, it is sometimes advantageous to inline a simple piece of code to tailor it to a specific application. While most iterative methods are too lengthy for inlining, conjugate gradients has a short implementation. The value of inlining pcg is that calls to afun and mfun may be inlined. Beware that in the following code, A must be symmetric positive definite—otherwise, it will not work!
%% Use PCG to tol = 1e−6; maxit = 20; x = x0; solve Ax = b %% % Convergence tolerance % Maximum number of iterations % Set x to initial guess (if no guess is available, use x0 = 0)
Compute w = A x
r = b − w; % r = residual vector
Set p = r, or if using a preconditioner, p = M −1 r
d = r'*p; for iter = 1:maxit if norm(p) < tol break; end
% PCG converged successfully (small change in x)
Compute w = A p
alpha = e/(p'*w); x = x + alpha*p; r = r − alpha*w;
Set w = r, or if using a preconditioner, w = M −1 r
if norm(w) < tol && norm(r) < tol break; % PCG converged successfully (small residual) end dlast = d; d = r'*w; p = w + (d/dlast)*p; end
After the loop, the final solution is in x. The method converged successfully with accuracy tol if it encountered a break statement, otherwise, it failed or needed more iterations. The relative residual may be computed by relres = norm(r)/norm(b) after the loop and resvec by storing the value of norm(r) each iteration.
21
9
Numerical Integration
b
Quadrature formulas are numerical approximation of integrals of the form f (x) dx ≈
a k
wk f (xk ),
where the xk are called the nodes or abscissas and wk are the associated weights. Simpson’s rule is
b
f (x) dx ≈
a
h [f (a) + 4f (a + h) + f (b)] , 3
h=
b−a . 2
h 4h h 3, 3 , 3.
Simpson’s rule is a quadrature formula with nodes a, a + h, b and node weights Matlab has several functions for quadrature in one dimension: quad quadl quadgk quadv adaptive adaptive adaptive adaptive Simpson Gauss-Lobatto Gauss-Kronrod Simpson
Better for low accuracy and nonsmooth integrands Higher accuracy with smoother integrands Oscillatory integrands and high accuracy Vectorized for multi-valued integrands
The quad- functions are robust and precise, however, they are not very efficient. They use an adaptive refinement procedure, and while this reduces the number of function calls, they gain little from vectorization and incur significant overhead. If an application requires approximating an integral whose integrand can be efficiently vectorized, using nonadaptive quadrature may improve speed.
for n = −20:20 % Compute Fourier series coefficients c(n + 21) = quad(@(x)(exp(sin(x).ˆ6).*exp(−1i*x*n)),0,pi,1e−4); end
This code runs in 5.16 seconds. In place of quad, using Simpson’s composite rule with N = 199 nodes yields results with comparable accuracy and allows for vectorized computation. Since the integrals are all over the same interval, the nodes and weights only need to be constructed once.
N = 199; h = pi/(N−1); x = (0:h:pi).'; w = ones(1,N); w(2:2:N−1) = 4;
w(3:2:N−2) = 2;
w = w*h/3;
% Nodes % Weights
for n = −20:20 c(n + 21) = w * ( exp(sin(x).ˆ6).*exp(−1i*x*n) ); end
This version of the code runs in 0.02 seconds (200 times faster). The quadrature is performed by the dot product multiplication with w. It can be further optimized by replacing the for loop with one vector-matrix multiply:
[n,x] = meshgrid(−20:20, 0:h:pi); c = w * ( exp(sin(x).ˆ6).*exp(−1i*x.*n) );
22
For this example, quadv can be used on a multi-valued integrand with similar accuracy and speed,
n = −20:20; c = quadv(@(x)exp(sin(x).ˆ6).*exp(−1i.*x.*n),0,pi,1e−4);
9.1 One-Dimensional Integration
b a
f (x) dx is approximated by composite Simpson’s rule with
= = = = (b − a)/(N−1); (a:h:b).'; ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; w * f(x); % Approximately evaluate the integral
h x w I
where N is an odd integer specifying the number of nodes. A good higher-order choice is 4th-order Gauss-Lobatto [4] (as used by quadl), based on
1 −1 1 1 1 f (x) dx ≈ 6 f (−1) + 5 f (− √5 ) + 5 f ( √5 ) + 1 f (1). 6 6 6
N h d x w I
= = = = = =
max(3*round((N−1)/3),3) + 1; % Adjust N to the closest valid choice (b − a)/(N−1); (3/sqrt(5) − 1)*h/2; (a:h:b).'; x(2:3:N−2) = x(2:3:N−2) − d; x(3:3:N−1) = x(3:3:N−1) + d; ones(1,N); w(4:3:N−3) = 2; w([2:3:N−2,3:3:N−1]) = 5; w = w*h/4; w * f(x); % Approximately evaluate the integral
The number of nodes N must be such that (N − 1)/3 is an integer. If not, the first line adjusts N to the closest valid choice. It is usually more accurate than Simpson’s rule when f has six continuous derivatives, f ∈ C 6 (a, b). A disadvantage of this nonadaptive approach is that the accuracy of the result is only indirectly controlled by the parameter N. To guarantee a desired accuracy, either use a generously large value for N or, if possible, determine the error bounds [6, 7] Simpson’s rule error 4th-order Gauss-Lobatto error ≤ ≤ (b − a)h4 max f (4) (η) provided f ∈ C 4 (a, b), a≤η≤b 180 27(b − a)h6 max f (6) (η) provided f ∈ C 6 (a, b), a≤η≤b 56000
b−a where h = N −1 . Note that these bounds are valid only when the integrand is sufficiently differentiable: f must have four continuous derivatives for the Simpson’s rule error bound, and six continuous derivatives for 4th-order Gauss-Lobatto.
23
Composite Simpson’s rule is a fast and effective default choice. But depending on the situation, other methods may be faster and more accurate:
If the integrand is expensive to evaluate, or if high accuracy is required and the error bounds above are too difficult to compute, use one of Matlab’s adaptive methods. Otherwise, consider composite methods. Use higher-order methods (like Gauss-Lobatto/quadl) for very smooth integrands and lower-order methods for less smooth integrands. Use the substitution u =
1 1−x
or Gauss-Laguerre quadrature for infinite-domain integrals like
∞ . 0
9.2 Multidimensional Integration
An approach for evaluating double integrals of the form a c f (x, y) dy dx is to apply one-dimensional b quadrature to the outer integral a F (x) dx and then for each x use one-dimensional quadrature over d the inner dimension to approximate F (x) = c f (x, y) dy. The following code does this with composite Simpson’s rule with Nx×Ny nodes:
%%% Construct Simpson nodes and weights over x %%% h = (b − a)/(Nx−1); x = (a:h:b).'; wx = ones(1,Nx); wx(2:2:Nx−1) = 4; wx(3:2:Nx−2) = 2; %%% Construct Simpson nodes and weights over y %%% h = (d − c)/(Ny−1); y = (c:h:d).'; wy = ones(1,Ny); wy(2:2:Ny−1) = 4; wy(3:2:Ny−2) = 2; %%% Combine for two−dimensional integration %%% [x,y] = meshgrid(x,y); x = x(:); y = y(:); w = wy.'*wx; w = w(:).'; I = w * f(x,y); % Approximately evaluate the integral
b d
wx = w*h/3;
wy = w*h/3;
Similarly for three-dimensional integrals, the weights are combined with
[x,y,z] = meshgrid(x,y,z); w = wy.'*wx; w = w(:)*wz; x = x(:); y = y(:); w = w(:).'; z = z(:);
When the integration region is complicated or of high dimension, Monte Carlo integration techniques are appropriate. The disadvantage is that an N -point Monte Carlo quadrature has error on the order O( √1 ), so many points are necessary even for moderate accuracy. Suppose that N points, N x1 , x2 , . . . , xN , are uniformly randomly selected in a multidimensional region (or volume) Ω. Then f dV ≈
Ω Ω
dV N
N
f (xn ).
n=1
24
To integrate a complicated region W that is difficult to sample uniformly, find an easier region Ω that contains W and can be sampled [5]. Then f dV =
W Ω
f · χW dV ≈
Ω
dV N
N
f (xn )χW (xn ),
n=1
χW (x) =
1, 0,
x ∈ W, x ∈ W. /
χW (x) is the indicator function of W : χW (x) = 1 when x is within region W and χW (x) = 0 otherwise. Multiplying the integrand by χW sets contributions from outside of W to zero. For example, consider finding the center of mass of the shape W defined by cos 2
2 2
x2 + y 2 x ≤ y
and x + y ≤ 4. Given the integrals M = W dA, Mx = W x dA, and My = W y dA, the center of M mass is ( Mx , My ). The region is contained in the rectangle Ω defined by −2 ≤ x ≤ 2 and −2 ≤ y ≤ 2. M The following code estimates M , Mx , and My :
%%% Uniformly randomly sample points (x,y) in Ω %%% x = 4*rand(N,1) − 2; y = 4*rand(N,1) − 2; %%% Restrict the points to region W %%% i = (cos(2*sqrt(x.ˆ2 + y.ˆ2)).*x <= y & x.ˆ2 + y.ˆ2 <= 4); x = x(i); y = y(i); %%% Approximately evaluate the integrals %%% area = 4*4; % The area of rectangle Ω M = (area/N) * length(x); Mx = (area/N) * sum(x); My = (area/N) * sum(y);
2
0
-2 -2
0
2
Region W sampled with N = 1500. The center of mass is ≈ (0.5, 0.7).
More generally, if W is a two-dimensional region contained in the rectangle defined by a ≤ x ≤ b and c ≤ y ≤ d, the following code approximates W f dA:
x y i x = = = = a + (b−a)*rand(N,1); c + (d−c)*rand(N,1); logical(indicatorW(x,y)); x(i); y = y(i);
area = (b−a)*(d−c); I = (area/N) * sum(f(x,y));
% Approximately evaluate the integral
where indicatorW(x,y) is the indicator function χW (x, y) for region W . For refinements and variations of Monte Carlo integration, see for example [1].
25
10
Signal Processing
Even without the Signal Processing Toolbox, Matlab is quite capable in signal processing computations. This section lists code snippets to perform several common operations efficiently.
Moving-average filter
To compute an N-sample moving average of x with zero padding:
y = filter(ones(N,1)/N,1,x);
For large N, it is faster to use
y = cumsum(x)/N; y(N+1:end) = y(N+1:end) − y(1:end−N);
Locating zero-crossings and extrema
To obtain the indices where signal x crosses zero:
i = find(diff(sign(x))); % The kth zero−crossing lies between x(i(k)) and x(i(k)+1)
Linear interpolation can be used for subsample estimates of zero-crossing locations:
i = find(diff(sign(x))); i = i − x(i)./(x(i+1) − x(i)); % Linear interpolation
Since local maximum and minimum points of a signal have zero derivative, their locations can be estimated from the zero-crossings of diff(x), provided the signal is sampled with sufficiently fine resolution. For a coarsely sampled signal, a better estimate is
iMax = find(sign(x(2:end−1)−x(1:end−2)) + sign(x(2:end−1)−x(3:end)) > iMin = find(sign(x(2:end−1)−x(1:end−2)) + sign(x(2:end−1)−x(3:end)) < ... 0) + 1; ... 0) + 1;
FFT-based convolution
This line performs FFT-based circular convolution, equivalent to y = filter(b,1,x) except near the boundaries, provided that length(b) < length(x):
y = ifft(fft(b,length(x)).*fft(x));
For FFT-based zero-padded convolution, equivalent to y = filter(b,1,x),
N = length(x)+length(b)−1; y = ifft(fft(b,N).*fft(x,N)); y = y(1:length(x));
In both code snippets above, y will be complex-valued, even if x and b are both real. To force y to be real, follow the computation with y = real(y). If you have the Signal Processing Toolbox, it is faster to use fftfilt for FFT-based, zero-padded filtering. 26
Noncausal filtering and other boundary extensions
For its intended use, the filter command is limited to causal filters, that is, filters that do not involve “future” values to the right of the current sample. Furthermore, filter is limited to zero-padded boundary extension; data beyond the array boundaries are assumed zero as if the data were padded with zeros. For two-tap filters, noncausal filtering and other boundary extensions are possible through filter’s fourth initial condition argument. Given a boundary extension value padLeft for x(0), the filter yn = b1 xn + b2 xn−1 may be implemented as
y = filter(b,1,x,padLeft*b(2));
Similarly, given a boundary extension value padRight for x(end+1), the filter yn = b1 xn+1 + b2 xn is implemented as
y = filter(b,1,[x(2:end),padRight],x(1)*b(2));
Choices for padLeft and padRight for various boundary extensions are Boundary extension Periodic Whole-sample symmetric Half-sample symmetric Antisymmetric padLeft x(end) x(2) x(1) 2*x(1)-x(2) padRight x(1) x(end-1) x(end) 2*x(end)-x(end-1)
It is in principle possible to use a similar approach for longer filters, but ironically, computing the initial condition itself requires filtering. To implement noncausal filtering and filtering with other boundary handling, it is usually fastest to pad the input signal, apply filter, and then truncate the result.
Upsampling and Downsampling
Upsample x (insert zeros) by factor p:
y = zeros(length(x)*p−p+1,1); y(1:p:length(x)*p) = x; % For trailing zeros, use y = zeros(length(x)*p,1);
Downsample x by factor p, where 1 ≤ q ≤ p:
y = x(q:p:end);
27
Haar Wavelet
This code performs K stages of the Haar wavelet transform on the input data x to produce wavelet coefficients y. The input array x must have length divisible by 2K . Forward transform:
y = x(:); N = length(y); for k = 1:K tmp = y(1:2:N) + y(2:2:N); y([1:N/2,N/2+1:N]) = ... [tmp;y(1:2:N) − 0.5*tmp]/sqrt(2); N = N/2; end
Inverse transform:
x = y(:); N = length(x)*pow2(−K); for k = 1:K N = N*2; tmp = x(N/2+1:N) + 0.5*x(1:N/2); x([1:2:N,2:2:N]) = ... [tmp;x(1:N/2) − tmp]*sqrt(2); end
Daubechies 4-Tap Wavelet
This code implements the Daubechies 4-tap wavelet in lifting scheme form [2]. The input array x must have length divisible by 2K . Filtering is done with symmetric boundary handling. Forward transform:
U = 0.25*[2−sqrt(3),−sqrt(3)]; ScaleS = (sqrt(3) − 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2); y = x(:); N = length(y); for k = 1:K N = N/2; y1 = y(1:2:2*N); y2 = y(2:2:2*N) + sqrt(3)*y1; y1 = y1 + filter(U,1,... y2([2:N,max(N−1,1)]),y2(1)*U(2)); y(1:2*N) = [ScaleS*... (y2 − y1([min(2,N),1:N−1]));ScaleD*y1]; end
Inverse transform:
U = 0.25*[2−sqrt(3),−sqrt(3)]; ScaleS = (sqrt(3) − 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2); x = y(:); N = length(x)*pow2(−K); for k y1 y2 y1 = 1:K = x(N+1:2*N)/ScaleD; = x(1:N)/ScaleS + y1([min(2,N),1:N−1]); = y1 − filter(U,1,... y2([2:N,max(N−1,1)]),y2(1)*U(2)); x([1:2:2*N,2:2:2*N]) = [y1;y2 − sqrt(3)*y1]; N = 2*N; end
To use periodic boundary handling rather than symmetric boundary handling, change appearances of [2:N,max(N-1,1)] to [2:N,1] and [min(2,N),1:N-1] to [N,1:N-1].
28
11
Miscellaneous
Clip a value without using if statements
To clip (or clamp, bound, or saturate) a value to within a range, the straightforward implementation is
if x < x = elseif x = end lowerBound lowerBound; x > upperBound upperBound;
Unfortunately, this is slow. A faster method is to use the min and max functions:
x = max(x,lowerBound); x = min(x,upperBound); % Clip elements from below, force x >= lowerBound % Clip elements from above, force x <= upperBound
Moreover, it works per-element if x a matrix of any size.
Convert any array into a column vector
It is often useful to force an array to be a column vector, for example, when writing a function expecting a column vector as an input. This simple trick will convert the input array of any shape and number of dimensions to a column vector, even if the input is already a column vector.
x = x(:); % Convert x to a column vector
Similarly, x = x(:,:) reduces an N-d array to 2D, where the number of rows stays the same.
Find the min/max of a matrix or N-d array
Given a matrix input (with no other inputs), the min and max functions operate along the columns, finding the extreme element in each column. To find the extreme element over the entire matrix, one way for a two-dimensional matrix is min(min(A)). A better method that works regardless of the number of dimensions and determines the extreme element’s location is
% Find the minimum element in A % The minimum value is MinValue, the index is MinIndex MinSub = ind2sub(size(A), MinIndex); % Convert MinIndex to subscripts [MinValue, MinIndex] = min(A(:));
The minimum element is A(MinIndex) or A(MinSub(1), MinSub(2), . . .) as a subscript reference. (Similarly, replace min with max for maximum value.)
29
Flood filling
Flood filling, like the “bucket” tool in image editors, can be elegantly written as a recursive function:
function I = flood1(I,c,x,y) % Flood fills image I from point (x,y) with color c. c2 = I(y,x); I(y,x) = c; if if if if x x y y > < > < 1 size(I,2) 1 size(I,1) & & & & I(y,x−1) I(y,x+1) I(y−1,x) I(y+1,x) == == == == c2 c2 c2 c2 & & & & I(y,x−1) I(y,x+1) I(y−1,x) I(y+1,x) ˜= ˜= ˜= ˜= c, c, c, c, I I I I = = = = flood1(I,c,x−1,y); flood1(I,c,x+1,y); flood1(I,c,x,y−1); flood1(I,c,x,y+1); end end end end
Being a highly recursive function, this is inefficient in Matlab. The following code is faster:
function I = flood2(I,c,x,y) % Flood fills image I from point (x,y) with color c. LastFlood = zeros(size(I)); Flood = LastFlood; Flood(y,x) = 1; Mask = (I == I(y,x)); FloodFilter = [0,1,0; 1,1,1; 0,1,0]; while any(LastFlood(:) ˜= Flood(:)) LastFlood = Flood; Flood = conv2(Flood,FloodFilter,'same') & Mask; end I(find(Flood)) = c;
The key is the conv2 two-dimensional convolution function. Flood filling a 40×40-pixel region takes 1.168 seconds with flood1 and 0.067 seconds with flood2.
Vectorized use of set on GUI objects
While not a speed improvement, a trick for cleaner code is vectorized use of set on GUI objects. A serious graphical user interface (GUI) can have dozens of objects and many properties to configure. For example, to define three edit boxes with white text background and left text alignment:
uicontrol('Units', 'normalized', 'Position', [0.1,0.9,0.7,0.05], ... 'HorizontalAlignment', 'left', 'Style', 'edit', 'BackgroundColor', [1,1,1]); uicontrol('Units', 'normalized', 'Position', [0.1,0.8,0.7,0.05], ... 'HorizontalAlignment', 'left', 'Style', 'edit', 'BackgroundColor', [1,1,1]); uicontrol('Units', 'normalized', 'Position', [0.1,0.7,0.7,0.05], ... 'HorizontalAlignment', 'left', 'Style', 'edit', 'BackgroundColor', [1,1,1]);
30
A vectorized call to set reduces this to
h(1) = uicontrol('Units', 'normalized', 'Position', [0.1,0.9,0.7,0.05]); h(2) = uicontrol('Units', 'normalized', 'Position', [0.1,0.8,0.7,0.05]); h(3) = uicontrol('Units', 'normalized', 'Position', [0.1,0.7,0.7,0.05]); set(h, 'HorizontalAlignment', 'left', 'Style', 'edit','BackgroundColor', [1,1,1]);
12
Further Reading
For more reading on vectorization, see the MathWorks vectorization guide: http://www.mathworks.com/support/tech-notes/1100/1109.shtml For further details and examples on JIT Acceleration, see for example http://www.mathworks.com/access/helpdesk_r13/help/techdoc/matlab_prog/ch7_perf.html For a thorough guide on efficient array manipulation, see MATLAB array manipulation tips and tricks: http://home.online.no/~pjacklam/matlab/doc/mtt For numerical methods with Matlab, see Numerical Computing with MATLAB : http://www.mathworks.com/moler
MEX
In a coding situation that cannot be optimized any further, keep in mind that the Matlab language is intended primarily for easy prototyping rather than high-speed computation. In some cases, an appropriate solution is Matlab Executable (MEX) external interface functions. With a C or Fortran compiler, it is possible to produce a MEX function that can be called from within Matlab in the same as an M-function. The typical speed improvement over equivalent M-code is easily ten-fold. Installations of Matlab include an example of MEX under directory <MATLAB>/extern/examples/mex. In this directory, a function “yprime” is written as M-code (yprime.m) and equivalent C (yprime.c) and Fortran (yprime.f) MEX source code. For more information, see the MathWorks MEX-files guide http://www.mathworks.com/support/tech-notes/1600/1605.html For information on compiling MEX files with the GNU compiler collection (GCC), see http://gnumex.sourceforge.net Beware, MEX is an ambitious alternative to M-code. Writing MEX files requires learning some low-level details about Matlab, and takes more time to develop than M-code.
31
Matlab Compiler
The Matlab Compiler generates a standalone executable from M-code. For product information and documentation, see http://www.mathworks.com/access/helpdesk/help/toolbox/compiler
References
[1] A. Bielajew. “The Fundamentals of the Monte Carlo Method for Neutral and Charged Particle Transport.” University of Michigan, class notes. Available online at http://www-personal.engin.umich.edu/~bielajew/MCBook/book.pdf [2] I. Daubechies and W. Sweldens. “Factoring Wavelet Transforms into Lifting Steps.” Sep. 1996. [3] J. W. Demmel. Applied Numerical Linear Algebra. SIAM, 1997. [4] W. Gander and W. Gautschi. “Adaptive Quadrature–Revisited,” BIT, vol. 40, pp. 84-101, 2000. [5] W. Press, B. Flannery, S. Teukolsky and W. Vetterling. Numerical Recipes. Cambridge University Press, 1986. [6] E. Weisstein. “Lobatto Quadrature.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/LobattoQuadrature.html [7] E. Weisstein. “Simpson’s Rule.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/SimpsonsRule.html
32
