Effective Programming and Data Types in MatlabCenter for Interdisciplinary Research and ConsultingDepartment of Mathematics and StatisticsUniversity of Maryland, Baltimore Countywww.umbc.edu/circWinter 2008Mission and Goals: The Center for Interdisciplinary Research and Consulting(CIRC) is a consulting service on mathematics and statistics provided by the Depart-ment of Mathematics and Statistics at UMBC. Established in 2003, CIRC is dedicatedto support interdisciplinary research for the UMBC campus community and the pub-lic at large. We provide a full range of consulting services from free initial consultingto long term support for research programs.CIRC offers mathematical and statistical expertise in broad areas of applications,including biological sciences, engineering, and the social sciences. On the mathematicsside, particular strengths include techniques of parallel computing and assistance withsoftware packages such as MATLAB and COMSOL Multiphysics (formerly known asFEMLAB). On the statistics side, areas of particular strength include Toxicology,Industrial Hygiene, Bioequivalence, Biomechanical Engineering, Environmental Sci-ence, Finance, Information Theory, and packages such as SAS, SPSS, and S-Plus.Copyrightc 2003–2008 by the Center for Interdisciplinary Research and Consult-ing, Department of Mathematics and Statistics, University of Maryland, BaltimoreCounty. All Rights Reserved.This tutorial is provided as a service of CIRC to the community for personal usesonly. Any use beyond this is only acceptable with prior permission from CIRC.This document is under constant development and will periodically be updated. Stan-dard disclaimers apply.Acknowledgements: We gratefully acknowledge the efforts of the CIRC researchassistants and students in Math/Stat 750 Introduction to Interdisciplinary Consultingin developing this tutorial series.MATLAB is a registered trademark of The MathWorks, Inc., www.mathworks.com.31 Optimizing M-files: VectorizationIn this section, we discuss vectorization – the main technique used in optimizingMatlab code. In general, vectorization entails carrying out computations on datastored as vectors (or matrices) using either linear algebra capabilities of Matlab orMatlab’s built-in functions that operate on vectors. As a general rule, loops are slowand one should instead utilize Matlab’s extensive matrix/vector capabilities wheneverpossible; the reason is that Matlab’s matrix/vector operations are fully optimized.Moreover, readers can be refereed to MathWorks’ website for a more detailed guidefor code vectorization. Its web page is onhttp://www.mathworks.com/support/tech-notes/1100/1109.htmlWe have already seen examples of vectorization in our discussion of logical sub-scripting. For example, the problem in which we wanted to zero out the elements ofa given matrix which fell below a given tolerance could have also been programmedusing a double for loop, which would also work but would be slower; the use of logicalsubscripting in that example was an example vectorization. Here we discuss furtheroptions to vectorize Matlab code.Example 1. Our first example of vectorizing a piece of code uses again the ideaof logical subscripting. Say we are given two n × n matrices A and B and we wantto form a matrix C which is defined by Cij= max(Aij, Bij); one way to solve thisproblem is to proceed as follows:C = zeros(n);for i = 1 : nfor j = 1 : nC(i,j) = max(A(i,j), B(i,j));endendHowever, the following (vectorized) version is both shorter and more efficient:C = B;I = A>B;C(I)=A(I);Example 2. Another example that shows exactly the same idea of logic sub-scripting is as follow. Give n a 200 × 400 matrix, we want to find out all entries thatare smaller than zero, and set them to be zero. You can prepare such matrix A bytyping A=rand(200,400)-0.5;. The non-vectorized code takes much time to finishthe problem:tic[m,n]=size(A)for i=1:m4 1 OPTIMIZING M-FILES: VECTORIZATIONfor j=1:nif A(i,j)<0A(i,j)=0;endendendtocThis is what C language programming does. The vectorized code does it conciselyand efficiently.ticA(A<0)=0;tocThe elapsed time are 0.1600 and 0.0500, respectively (on the same machine). Isn’t itshocking? The above examples show the spirit of vectorization in general; we wouldlike to replace loops by operations that utilize Matlab’s matrix/vector capabilities.The next example is more mathematical in nature.Example 3. Given two vectors a and b, one defines their tensor product a ⊗ b(a matrix) by the following(a ⊗ b)ij= aibj.The following Matlab function returns the tensor product of the vectors a and b:function C=tensor1(a,b)n = length(a);C = zeros(n);for i = 1 : nfor j = 1 : nC(i,j) = a(i)*b(j);endendThe following vectorized code does the same thing much more efficiently (at thesame time a much shorter code too):function C=tensor2(a,b)C = a(:)*b(:).’; % column * rowHere, a(:) converts any row or column vector into a column. b(:).’ is the non-conjugate transpose. For more information, type help transpose, and help ctranspose.You can test these two codes like below.>>a=[1:1000];>>b=[1000:-1:1];>>tic; tensor1(a,b); toc>>tic; tensor2(a,b); toc5Function Descriptionmin Find the smallest componentmax Find the largest componentsum Find the sum of array elementscumsum Find cumulative sumfind Find indices and values of nonzero elementsall Test to determine if all elements are nonzeroany Test for any non-zerosprod Find product of array elementscumprod Find the cumulative product of array elementsrepmat Replicate and tile an arrayreshape Change the shape of an arraysort Sort array elements in ascending or descending orderunique Find unique elements of a setTable 1: Some of the Matlab’s built-in functions used in vectorizationThe elapsed time are 0.1300 and 0.0300, respectively (on the same machine). Moresophisticated vectorization uses Matlab’s built-in function which operate on vectors.Table ?? lists some of the most commonly used Matlab functions in vectorization.Example 4. Consider the operation of computing the scalar product of two n×nmatrices with real entries. We define,A · B =nXi=1nXj=1AijBij.The Matlab function testmatrixdot in the next page provides two implemen-tations of the the above operation in the subfunction scalar_for and scalar_vec.The function testmatrixdot receives n as an input parameter, generates two n × nrandom matrices A and B and computes the wall clock run time of both vectorizedand non-vectorized ve rsion for